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Interferometric test for -wave flow 
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Ballistic Research Laboratories, Aberdeen Proving Ground, \Marvland 


(Received 8 December 1956) 


SUMMARY 

An analytical expression for the fringe shift in the V-wave 
is derived from the improved linearized theory for a slender 
supersonic projectile. From this expression an approximate 
mapping function is found which gives a simple test for N-wave 
flow. The validity of the fringe-shift expression is quantitatively 
confirmed by measuring an interferogram of the flow around a 
sphere. N-wave flow is shown to exist around a small sphere 
for r greater than about 70 diameters. Measurements of the 
shock waves from a sphere and a cone-cylinder show that the 
shocks assume their asymptotic positions for r > 14 diameters 
for the sphere, and r > 7 diameters for the cone-cylinder. 


1. INTRODUCTION 
1.1. Interferometric methed of analysis 


Interferometric investigation of the fluid flow around a supersonic 
projectile in free flight gives a quantitative record of the density over the 
entire flow field (Ladenburg & Bershader 1954). Experience shows that 
the reduction of the fringe shift to density values for an axisymmetric flow 
is a cumbersome and time-consuming process because, in general, the 
relationship between fringe shift and density at a particular point is not 
simple. It is useful, therefore, to find flow regions for which some intrinsic 
property can be determined directly from the measurements of fringe 
shift. 

An interferogram of a supersonic cone-cylinder shows two regions 
which have conspicuous symmetry of fringe shape, suggesting the possibility 
of a simplified analysis. ‘The first of these is the region near the cone in 
which the fringes are nearly straight and parallel (Giese, Bennett & Bergdolt 
1950; Giese & Bergdolt 1953; Bergdolt 1953; Cole, Solomon & Willmarth 
1953). ‘The flow here is characterized by the fact that the physical variables 
are constant along straight lines through the vertex. Assumption of this 
flow regularity leads to a method of plotting fringe shift which verifies in 
many instances the approximation of real flows to idealized conical flow. 

The second of these regions lies between the front and rear shock waves 
at rather large distances from the projectile axis. Here the fringes have a 
gentle curvature and a similarity of shape which change only slowly as the 
F.M. 
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distance from the axis increases. ‘This similarity suggests an underlying 
simplicity of the fundamental flow field. Experiments by DuMond et al. 
(1946) showed that at large radial distances, the pressure profile parallel 
to the axis consists of a sudden rise at the front shock, followed by a linear 
decrease to a value below that for the free stream, and then a sudden rise 
at the rear shock. The curve so generated has the shape of a capital N; 
hence is given the name ‘ .V-wave’. 


1.2. Scope of the paper 
Using the results obtained by Whitham (1952) in his improved linearized 


theory for slender supersonic projectiles, we derive here an analytical - 


expression for the fringe shift in the N-wave region. We investigate this 
function and find a mapping. property which gives a simple criterion for 
comparing theoretical and measured data. 

Interferograms taken of the flow around spheres and slender cone- 
cylinders are measured, and observed fringe shift is compared with the 
theory. ‘To justify the use of spheres, we discuss the validity of comparing 
the distant flow around blunt bodies with values predicted by slender-body 
theory. Figure 1 (plate 1) shows a general view of the flow about a sphere. 


2. "THEORY 
2.1. The fringe-shift integral 
If (x,r) are the cylindrical polar coordinates of the axisymmetric dis- 
turbance, and the front tip of the projectile is at the origin with the line 
of flight parallel to the x-axis, then the fringe shift 6(x, 7) is related (Bennett 
et al. 1952) to the density p at (x,r) by 
(p—po)t dt (2. 
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Figure 2. Path of integration of the fringe-shift integral for a supersonic projectile. 
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where K is the Gladstone—Dale constant, A the wavelength of the light in 
vacuum, py the free-stream density, ry the outer radius of the disturbance 
at x, and ¢ is the variable of integration in the r direction. As in figure 2, 
the fringe shift at (x, r) is determined by an integration of the density values 
found along the line from r to ry. 

An inversion of (2.1) allows the reduction of fringe-shift measurements 
to density throughout the flow field. Such a reduction does not interest us 
here; for our purpose is to obtain in the N-wave a functional description 
of (p — py) and from it derive an expression for the fringe shift in that region. 


2.2. Improved linearized theory 

The linearized theory of Karman & Moore (1932) for a slender supersonic 
projectile gives solutions which are good first approximations to the actual 
conditions at the surface of the projectile, but which fail as the distance r 
increases. ‘This failure arises from the fact that the Mach lines are parabolae 
in a second approximation, while the linear approximations to them are 
straight lines; and although a curved Mach line intersects a straight one 
at the surface of the body, the curves diverge with increasingr. Theimproved 
theory offered by Whitham (1952) makes the linear solutions uniformly 
valid over the entire field by associating them with the more exact Mach 
lines. ‘This improved theory is still based upon a first approximation to 
the potential flow equations, and, like linear theory, neglects terms of order . 
u*, v*, where uw and v are the small perturbation velocities in the x and r 
directions. Nevertheless, it has the additional advantage over linear theory 
that the existence and position of the shock waves are predicted. 

With this background in mind, we now discuss that part of the Whitham 
theory necessary to our development. Equations which appear in Whitham’s 
paper are identified by the notation Wh. 

1. The shock waves. At great distances from its axis, the supersonic 
projectile produces two shock waves, both extending to infinity. The 
equations for the shock waves at large r are, from Wh(43), 


where the upper sign represents the front shock and the lower sign the rear 
shock, x? = M?—1, and A and yp are constants related to body shape. 
The straight line « = ar + yp and the lines given by (2.2) intersect the x-axis 
at yy. By inspection of (2.2) one sees that the front shock wave lies ahead, 
and the rear shock behind the straight line by the amount Ar'*. Thus, 
the horizontal distance between the two shocks is 2A4r!*. This property 
will be used later to determine A experimentally. The fact that the locus 
of the mid-points of the horizontal distances is the line x = ar+ypo will 
enable us to determine « when the projectile axis is known, or the projectile 
axis when « is known. 

2. The pressure distribution. When r is large, the pressure distribution 
between the shock waves is by Wh(71), 


(P —Po)/Po = yM*(2a) — + Yo) (2.3). 
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where & = 2°-1*(y+1)M'x 32, p is the disturbed and p, the free-stream 
pressure. Notice that along a trace of constant r the pressure difference p — p, 
decreases linearly with « from a positive value at the front shock to an equal 
negative value at the rear shock. Midway between the shocks on the 
straight characteristic, we have p = py. The pressure slope, y.M?(kr)-1(2«)-?? 
depends only on the flow constants and the distance from the axis, and not 
upon the shape of the body producing the disturbance. In addition to 
this expression for the .V-wave distribution, which only exists in the distant 
regions, Whitham derives a more general equation which gives the pressure 
distribution at any distance from the projectile axis. Figure 3 shows a 
typical pressure signature along a horizontal trace for three distances from 
the axis :—near to the body and before the rear shock is formed, an intermediate 
position with both shocks, and the final wave. 


Figure 3. Typical pressure signatures along horizontal traces: top curve, near the 
body; middle curve, at an intermediate distance; bottom curve, in the 
N-wave region. 


2.3. Density distribution in the N-wave 

We now relate pressure and density in the .\-wave. Despite the seeming 
simplicity of this programme, some care must be exercised to assure an 
approximation consistent with that of linearized theory. 

It can be shown (Heaslet & Lomax 1954) that the entropy may be 
considered constant across a shock wave if third and higher-order powers 
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of the perturbation velocities are neglected. We neglect second and higher- 
powers and may, therefore, assume that the flow behind the shock is 
isentropic. ‘The adiabatic gas law applies; thus 


= (P/Po)*”: (2.4) 
Expanding (2.4) in Taylor series around p = py, we get 
—) =14+-(——)+ + . 
Y\ Po 2y? \ Po Po 
In 1Wh(66), Whitham finds that 
(P—Po) Po = Olu); 


so to the linear approximation 


Po Po Y\ Po 


Po Po 
If we substitute (2.3) into (2.5), we find that 
py = py MPR\(2x) — (2.6) 


This is the density distribution between the shock waves in the V-wave 
region. Along a trace of constant 7, the density profile of the N-wave has 
the same characteristic shape as the pressure curve. 


2.4. Fringe shift between the shock waves 
Combining (2.6) and (2.1), we obtain 


\2Kpy) M2 


Ty (af —x+y 
| dt. (2.7) 


Here x is constant along the path of integration from r to ry. Regrouping 
terms and integrating, we find 


Nk 


where ry is related to x by means of the front shock equation 


x= ory+yy— Art. (2.9) 


&(x,r) = 


7) = 


Equations (2.8) and (2.9) predict the fringe shift between the shock waves. 

Let us examine the variation of d(x,r) along a trace r = const. as x 
increases from xy to vg, which lie on the front-and rear shocks respectively. 
Also, No lies on the line x — zr = 9, halfway between x, and xz (see figure 4). 
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Let f(e) g(e) = logle+(e?—1)!?], and R =(x—yp)/ar, where 
€=/Fy/r. For a particular flow the magnitude and sign of 5 is determined 
from (2.8) by the term in the curly brackets, i.e. by {f— Rg}. 


Figure 4+. The region between the two shock waves. 


By comparison with the equation of the straight characteristic 


[<1] 
R\=1; if = 
NX. 


Thus R increases monotonically with x from a value less than 1 at the front 
shock to a value greater than 1 at the rear shock. From the equation of the 
front shock wave, 


< R <14(A/a)r4, 


The variation of f(e) and Rg(e) with « is shown in figure 5 for R < 1. 
‘The curves have vertical tangents at « = 1 and intersect at no other point. 
(A proof of this fact appears in Appendix I.) When R > 1, they are related 
as in figure 6. In addition to their intersection and common vertical tangent 
at «= 1, they intersect at one other point, say « = «, (see Appendix 1), 
where e, increases as R increases. 

We can now examine the combination of functions {f—Rg}. As 
increases, so do «, R, and e, (when R >1 ). 

When x, <x <x, R<1, so f>Rg. Therefore, 5 is positive in the 
forward portion of the N-wave. 

When x) <x R>1, and f= Rg ife When x is just greater 
than %, « ><«,, so 6 is positive. But it is possible that as x increases, 
€, increases faster than e. If it does, then the conditions e <e, can occur, 
and 6 will pass through zero to negative values near the rear shock wave. 
The proof that de,/dox > de/dx is too complicated to be attempted here. 
We resort to calculation to establish the behaviour of 6 in this region. 
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Figure 6. Variation of the functions f and Rg when R> 1. 
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A calculation of 6 vs x alongr = const. from (2.8), using representative 
values of A, yy and the flow constants, yields the curve in figure 7. The 
fringe shitt does go to negative values near the rear shock wave. 
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Figure 7. Variation of fringe shift along a horizontal trace, calculated from (2-8) 


2.5. Fringe-shift mapping law 

Let us investigate the behaviour of the 6 vs x curve of figure 7 at 
several radial distances. ‘lable 1 gives the calculated values of 6 along 
three horizontal traces. Preliminary calculations* of 6 vs x show that d 
becomes zero about three-quarters of the way between x, and vg, no matter 
what the value of 7. This indicates that the ratio (w—.x,.)/(vs—w%y) has a 
significance which is independent of r; so in table 1, 6 is given for values of 


= Definition 1 (2.10) 
rather than for values of x. ‘Theoretically, vg— x) = 2Ar'* by (2.2), and so 
= (w—xy)/2Ar"4, Definition (2.11) 


Between the shock waves, 0 <€ <1. ‘The variable € seems indeed to be 
independent of r in that all the curves of 6 vs € given in table 1 pass through 
zero and reach positive maxima for the same values of €. However, the 
magnitude of 6 decreases as r increases. 

Because of the inverse variation of 6 with r, it seemed possible that for 
any two traces the relation 6,7; = 6,73 might hold for each value of € along 
the traces. When values of 6 and r were substituted in this equation, 
n turned out to be very close to }. For each trace, then, the values of 5 
were multiplied by 7°. A plot of dr'° vs € yields a curve which remains 
closely the same for all horizontal traces through the N-wave. Comparison 


* A discussion of this investigation is being published in an Office of Ordnance 
Research Report on the Second Annual Conference of Ordnance Mathematicians, 
February 1956. Inquiries should be addressed to Office of Ordnance Research, 
Box CM, Duke Station, Durham, North Carolina. 
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- of the last three columns of table 1 shows the agreement possible. ‘The 
last two columns agree to better than 1°, for many of the entries. ‘The 
variables € and 6 are dimensionless; values of the coordinates are rendered 
dimensionless by dividing them by the projectile diameter. 


! 
8 
g 
r= 10 r = 50 r= 100 r= 10 r= 50 r = 100 
0 0 0 0 0 0 0 
0-1 0-840 2 0-676 0-618 1-120 4 1-102 1-099 
0-2 1-004 8 0-808 0-740 1-318 1-316 
0-3 1-009 0 0-812 0-742 1-345 5 1-324 1-320 
0-4 0-909 4 0-731 0-666 #2127 1-192 1-184 
0-5 0-732 0 0-585 0-534 0-976 1 0-954 0-950 
0-6 0-490 6 0-385 0-350 0-654 2 0-628 0-622 
0-7 0-196 6 0-145 0-130 0-262 2 0-237 0-231 
0:8 —0-1466 | —0-138 —0-130 —0-1955 | —0-225 —0-231 
0-9 —0-5136 | —0-455 —0-422 —0-7089 | —0-742 —0-750 
1-0 —0-9556 | —0-807 —0-742 —1-2743 | —1-316 —1-320 


Table 1. Fringe shift calculated from equation (2:8); flow constants are dA — 2, 
a = 2 Ve = 2: 


We now show that the 6r'* mapping law is contained implicitly in the 
exact N-wave fringe-shift equation (2.8), and can be revealed by making 
uitable approximations. 

Referring to § 2.4 and making use of the definitions given there, we can 
write (2.8) as 

= Cri f—Re}; (2.12) 
where C = (2x)!*Kp,) W?(Ak) 1, and f, g are functions only of € (= ry/r). 
Further, let (4 /«)r-4 = o, and observe that o is proportional to Whitham’s 
approximate expression for the shock strength, Wh(44). Our task is now 
to find R, f and g as functions of o and €. For R this is comparatively easy 
and leads to a general proof of the possibility for f and g. 

The geometric relationship between an external point (x,r) and the 
corresponding points on the front shock wave (x,,r) and (x,7ry) is made 
clear schematically in figure +. Selecting the minus sign in (2.2), we may 
write 

Xy = (2.13 a) 
and x= ory tyy— (2.13 b) 
Substituting (2.13 a) into (2.11), the definition of €, we obtain 
€ = (x—ar—yyt+ Art) /24r"4, 
and solving this for (v—4»9)/a gives 


R= (2.14) 
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Division of (2.13b) by xr yields a relation between ¢, o and R as follows: 
e—oel4-R= 0. (2.15) 


From the quartic equation (2.15), it is clear that « = e(o,R). From 
(2.14) it follows that ¢ is a function only of o and €; consequently, the 
same is true of f and g. It remains now to find suitable expressions for 
these two functions. 

As a preliminary step, we find bounds for «. From the definition it is 
clear that 1 <e <ry/ry. The shock-wave equations (2.2) evaluated at x 
yield 


ary — Arl4 = ar + Arlt, 
or upon dividing by arg, 
(ry/rs)—1 = + 1], (2.16) 


where o, = (A/x)rc**. As rx, the shock strength o, approaches zero, 
and it follows from (2.16) that r,/r,>1. By choosing r, sufficiently large, 
we can bring « as close to unity as we please. We may therefore set « = 1+», 
where v = v(ry) < 1 is as small as we please. 

The quantity v may now be evaluated by approximating the solution 
of the quartic (2.15). Saving terms up to O(v?) in the binomial expansion 
of «!4, we find 


v= (2.17) 
With (2.14) this becomes 
v = O(or?), (2.17 a) 


Equation (2.17 a) shows that v is of the order o€. By (2.10) or (2.11) it is 
clear that € = O(1); hence we conclude that v is O(c). Upon saving terms 
of order o” in (2.17 a), we have 


vy = 2of(1+ 40) + (2.18) 


By definition, f = (e?—1)!? = (2v+v*)#*. Using (2.18) and expanding 
the square root, we find that 


f = + Jot + (2.19) 


A suitable expression for g is slightly more difficult to get. We proceed 
as follows. Solve f = (e?—1)!? for «, and obtain « = (f?+1)'?. Then 


g = logle 1)'"] = (2.20) 


‘The expansion tor (2.20) may be obtained in the form 


g= (2.21 a) 
1.3.5...(2n—3) 
where b, =(-1) (2.21 b) 


g 
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and f? = (e?—1) <1. This convergence conditions requires that « < v2, 
which, from our discussion of (2.16) can always be fulfilled by making r 
large enough. 

Using (2.14), (2.19) and (2.21a,b) in (2.12), we find after some 
manipulation 


8 = — 48) + O(0°?). (2.22) 


Since o = (A/x)r34, the product ro®? yields a term in r'*, When 
transposed, this provides the 6r'* scaling law discovered empirically. 
Henceforth, we shall call the quantity r1°(2C)14-33? the ‘scale factor’, 


and denote it by the letter G. ‘Then (2.22) may be written* 
8G = £41 — 42), (2.23) 


an expression which not only embodies the r'* mapping law, but also 
exhibits the zero at € = 0-75 already expected from the early investigation. 


5G 
A=12 
A=2,0=2 A= hand 
h(é) 
0 0 0 0 0 0 0 0 
0-1 0-280 0-276 0-275 0-276 0-275 0-290 0-274 
0-2 0-335 0-329 0-329 0-330 0-348 0-328 
0-3 0-336 0-331 0-330 0-330 0-329 0-349 0-329 
0-4 0-303 0-298 0-296 0-298 0-299 0-314 0-295 
055 0-244 0-239 0-237 0-238 0-237 0-251 0-236 
0-6 0-164 0-157 0-156 0-157 0-156 0-168 0-155 
0-7 0-066 0-059 0-058 0-058 0-057 0-065 0-056 
0-8 — (0-049 0-056 0-058 —0-058 — 0-062 0-056 — 0-060 
0-9 0-177 0-186 —0-188 —0-187 | —0-190 0-192 | —0-190 
1-0 0-319 0-329 — 0-330 0-330 — 0-332 —0-340 0-333 


Table 2. Scaled fringe shift calculated from the exact form, equation (2.8), and 
from the approximate form, h(€) = €'? (1— 4); G = ‘scale factor *. 


We can get an idea of the amount of approximation involved in this 
equation by comparing values of 5G calculated from the exact fringe-shift 
expression (2.8) with values of A(é) = £2(1—é). The comparison is 
made in table 2 for several combinations of r, A, and x We notice that 
for constant A and «, the approximation improves as / increases; and 
for constant r, as the ratio A/a decreases. This is consistent with the fact 
that the approximate form A(€) is obtained by neglecting powers of the 
term o = (A/a)r-34, 


* The authors are indebted to Dr Raymond Sedney for helpful discussions of 
the approximation procedure. 
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For A/x << 1, the approximation is within 1% of the exact value for all 
traces with r > 50 diameters. For A/x = 1-7 andr = 50, the approximation 
is within 6°,, of the exact value. 

Equation (2.23) shows that 6G is constant along curves where 
€ = (x—x,)/2Ar!4 = const.; or since 

xy = Ar", 
along the lines 
x= art yo + Arl4(2é — 1). (2.24) 
‘Vhis family of curves, shown schematically as dotted lines in figure 5. 
includes the shock waves (when € =0 and 1) and the straight characteristic 
x= ar+yy (when = 0-5). 


x 


Figure 8. Lines of constant 6G. 


Equation (2.23) gives a convenient test for the existence of .\-wave flow 
no matter what the Mach number. The fringe shift 6 can be measured 
along a horizontal trace and multiplied by the scale factor. ‘The existence 
ot -wave flow is established if this quantity, plotted against €, lies closely 
along the curve A(€) = €'2(1 — 32). 


2.6. Application of theory to flow around blunt bodies 

Since all of the theory described above is based upon an V-wave pressure 
expression derived from slender-body theory, one would naturally select 
a slender projectile in any attempt to substantiate the theory experimentally. 
However, stability considerations, to be discussed later, make the slender 
projectile impractical and require that a sphere be used in the investigation. 

From a theoretical standpoint it seems possible to use the disturbance 
from a sphere to verify the predicted \-wave fringe shift. Whitham (1950) 
shows that the pressure profile between the shocks far enough away from 
any supersonic body is linear, and that the pressure equation and the 
equation of the front shock take the form 


— Po= +h—ar), (2.29) 
and x= ar—br!4—h, (2.30) 


respectively, where 6 and / are constants depending on body shape. ‘These 
equations are of the same form as (2.2) and (2.3), which resulted from 
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slender-body theory. Although 6 and h depend upon the body shape, 
they are not explicitly set out in terms of body shape; whereas the constants 
A and yp of (2.2) and (2.3) are. We infer, therefore, as Lighthill (1954) 
has done, that equations (2.2) and (2.3) from Whitham’s slender-body 
theory apply to the disturbance sufficiently far from the axis of any supersonic 
projectile, provided that the constants A and yp are evaluated empirically 
from the shock wave positions and other data. 


3. EXPERIMENTAL PROCEDURE 
3.1. Theory of the experiments 

‘The experimental investigation is designed (1) to test the validity of 
the fringe-shift equation (2.8) between the shock waves, (2) to determine 
how near to the projectile axis the conditions characteristic of an ideal 
\-wave may be observed, and (3) to determine how near the axis the actual 
shock waves lie along the lines described by (2.2). 

The first studies were made of the flow around a 0-225 in. cone-cylinder 
in a region about 15 projectile diameters from the axis (see first footnote 
of §2.5), The fringe shift was much greater than that predicted for the 
ideal N-wave, and, in order to observe the flow further from the axis, it 
became necessary to use projectiles of smaller diameter. Final results 
were obtained with a 1 in. diameter sphere, which permitted observations 
to be made at distances up to 70 diameters. 

To substantiate the fringe-shift equation, the actual fringe shift is 
measured along a horizontal trace, multiplied by the scale factor, and 
compared with the function A(é) = €'?(1— 36). 

The scale factor is a function of C, x, and r; however, since the 
disturbance is caused by a blunt body, A cannot be calculated from the 
body shape, but must be determined from measurements of the shock wave 
positions with the help of the definition 4 = (xyg—wxy)/2r’4. ‘The reliability 
of this method is tested by measuring the shock positions from a } in. 
diameter sphere. Average values of A and yy are calculated from these 
positions and are used with (2.2) to plot the theoretical N-wave shock 
positions. These are compared with the measured positions. 

From a theoretical point of view it would have been desirable to compare 
the measured shock positions of a projectile of known velocity and line of 
flight with the positions predicted by the shock equations (2.2). Actually 
it is impractical to determine the line of flight accurately in the free flight 
range. In order to fix the axis in an interferogram, the existence of the 
straight characteristic, w = «+9, lying halfway between the shocks, is 
assumed. The picture is adjusted until the slope of this line corresponds 
to the « calculated from velocity measurements. Our method for verifying 
the .V-wave shock equations is thus a semi-empirical one. 


3.2. Apparatus and instrumentation 

To obtain quantitative data the projectiles were fired from a calibre 
(:22 rifle approximately along the axis of the Controlled Temperature- 
Pressure Range of the Air Flow Branch, Exterior Ballistics Laboratory, 
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The flow was observed through optical glass windows at the interferometer 
station of the range by a Mach—Zehnder interferometer with an 8 in.x 10 in. 
working field. (For a detailed discussion of this type of instrument, see 
Ladenburg & Bershader 1954.) 

In order to record the disturbance as far from the axis of the projectile 
as possible, the rifle was aimed to place the projectile near the edge of the 
interferogram. Fringes were adjusted to lie parallel to the trajectory. 

The light source employed at the interferometer is an exploding wire 
with a duration of about 3 microseconds (Lewis & Sleator 1956). Since 
the projectiles travelled about 1-5 mm during this time, their images were 
stopped on the film with a rotating mirror (Bergdolt et a/. 1951). A band 
of light about 20 A wide centred at the 4358 A mercury line is isolated from 
the continuous spectrum of the exploding wire by a slit and a carbon 
disulphide liquid-filled prism (Sleator et al. 1952). 

The projectile velocity is measured between a pair of stations 16-89 ft. 
apart served by 0-1 megacycle chronograph counters triggered by impulses 
from photo-cells which respond to fluctuations in light screens through 
which the missile passes. ‘The impulse from the down-range station also 
triggers a delay circuit, which fires the exploding wire of the interferometer 
station at the proper time. 

The interferograms are measured with a ‘Telereader—lelecordex 
combination. The Telereader projects a magnified image of the interfero- 
gram on a screen; horizontal and vertical cross hairs can be moved 
independently across the screen to measure distances in the r and x directions 
respectively. ‘The distances moved by the cross hairs are displayed by the 
counters of the Telecordex and can be permanently recorded with an 
automatic typewriter. If the counters are zeroed when the cross hairs are 
at one point of the interferogram (for example, at the nose of the projectile), 
the coordinates of any other point are measured by moving the cross hairs 
to that point. The measurement is accurate to within two Telereader units. 


3.3. The slender cone-cylinder 

A cylinder of 0-225 in. diameter with a conical tip of 20° included angle 
was initially selected for this investigation. This projectile is stable in 
flight when fired with small initial yaw and high spin (cf. McShane et al. 
1953), and it is just slender enough to meet the requirements of slender-body 
theory. A comparison of the results of Whitham’s theory and the exact 
adiabatic theory of flow past a cone shows reasonably good agreement for 
included cone angles up to 20°. 

The cylindrical portion of the projectile was made slightly oversize and 
was used as the rotating band in order to reduce the shock waves which attach 
themselves to the forward edge of a conventional rotating band. Even so, 
compression waves are propagated from the grooves engraved on the 
cylinder. 

An interferogram (figure 9, plate 2) was taken of the flow around the 
projectile at a Mach number of 2-25. Fringe-shift measurements taken at 
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distances from the axis up to about 13 diameters indicated that the ideal 
\-wave had not formed within the field of view. The results of these 
measurements are reported in a later section only as a qualitative comparison 
with results obtained with spheres. 

To observe the flow still farther away from the projectile axis, un- 
successful attempts were made to launch 0-1 in. diameter cone-cylinders. 
Sub-calibre projectiles such as these are made to fit the rifle bore by encase- 
ment in a plastic holder, called a sabot, which separates from the projectile 
after leaving the muzzle of the rifle. Difficulties arise because during 
separation the sabot imparts a certain amount of angular momentum and 
initial yaw to the projectile, and causes a body of borderline stability (such 
as the 20° cone-cylinder) to yaw violently or tumble in flight. It was therefore 
necessary to use a more stable projectile than the cone-cylinder. The 
sphere is a practical choice. It is stable, easily saboted, and can be readily 
obtained in almost any size. 


3.4. The spheres 


The spheres used were standard steel balls with diameters of } in. and 
7,in. Since these were sub-calibre projectiles, sabots were necessary. 
The sabots were plastic cylinders 0-4 in. long and 0-225 in. in diameter 
whose forward faces were hollowed as inverted cones; the spheres were 
placed in shallow holes drilled at the apex of the cones. 

An interferogram was taken of the flow around a } in. diameter sphere 
(figure 10, plate 3) at a Mach number of 1-35. Shock-wave measurements 
from this interferogram are compared with the ideal N-wave shock equations 
(2.2) in § 4.2. Fringe-shift measurements, made at distances from the axis 
up to 40 diameters, indicated the need for a still smaller projectile. 

Final firings were conducted with the 4 in. sphere. The small size 
of this projectile presented two technical difficulties. (1) The muzzle 
velocity was somewhat unpredictable with the result that the delay set into 
the station circuits before a round was fired was not always compatible 
with the velocity attained. (2) The intensity of the light falling on the 
photocells at the velocity stations was frequently not influenced enough 
by the passing projectile to trigger the velocity counters. It was necessary, 
therefore, to fire many rounds to get two usable interferograms. For one 
of these interferograms the velocity of the round was measured, and the 
shock positions are clearly defined, but the quality of the fringes is too poor 
to allow accurate fringe-shift measurement. On the other hand, the fringes 
in the other are distinct and measurable, but the velocity was not obtained. 
No change was made in the inclination of the fringes nor in the gun position 
between the firing of these two rounds. Thus the inclination of the axis 
of flight with the fringes is presumably the same in both interferograms. 
Fringe-shift measurements are made on the interferogram shown in 
figure 11 (plate 4). The projectile axis is determined in conjunction with 
the interferogram for which the velocity is known. | 
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+. INTrERFEROGRAM MEASUREMENT AND CALCULATIONS 


4.1. Summary of the data 

‘Table 3 gives a summary of measured and calculated data. Items 
1 to 8 list data measured in the laboratory at the time of firing. The air 
density, item 9, is determined from appropriate handbook tables using the 
measured values of temperature, pressure and humidity. The projectile 
velocity and Mach number, items 10 and 11, are calculated quantities; a 
discussion of velocity calculation for spheres is given in Appendix II. 
The quantity « is calculated from the Mach number for the projectile of 
figure 10(plate 3); for the + in. sphere of figure 11 (plate 4), « is the measured 
slope of the straight characteristic. ‘The constants 4 and yp, items 13 and 
14, are calculated from shock-position measurements in a manner described 
in the next paragraph. ‘To determine the undisturbed fringe spacing , 
the distances between adjacent fringes are measured in the undisturbed 
region in front of the shock wave; A is an average of these measurements. 

For those quantities which are used in a quantitative comparison with 
theory, table 3 gives the estimated error in the quantity. Where the datum 
is the average of a number of measurements or calculations (items 13, 1+ 
and 15), the estimated error is given in the form of a standard deviation. 


4.2. Verification of the shock positions 

Before any measurements of the interferogram can be made, the picture 
must be oriented on the Telereader so that the origin of the coordinate 
system lies at the nose of the projectile and the x-axis along the axis of 
Hight. Although this axis is not known, the velocity (and therefore ~) 
has been measured. If the slope of the straight line running midway 
between the shocks is made to have the value « by adjusting the picture 
of the Telereader, then the x-axis of the Telereader will lie along the 
projectile axis. 

The interferogram is initially placed in an arbitrary position on the 
‘Telereader. For some value of r, xy) and x, are measured and x, calculated 
from xy = (x,+xy); for another value of 7, another x) is determined in 
the same way. Now the coordinates of two points are known, and both 
points should lie on the line « = ar+ yo if the interferogram is properly 
oriented. Therefore, the slope of the line between the two points is 
calculated, and compared with the value of « calculated from the known 
velocity. If there is a discrepancy, the picture is rotated on the Telereader, 
and the same procedure is repeated. When the slope of the straight 
characteristic and x are as near to the same value as the Telereader 
adjustment will allow, then the x-axis is presumed to lie along the projectile 
axis. 

Now the coordinates of the shock waves (x\,7) and (x,,7) are measured 
for many values of r. For each value of r we calculate (1) the point (x), ’) 
on the straight characteristic from the relation x) = $(wy+.x,); (2) yp from 
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Figure 1. Single-fringe interferogram of a 0-218 in. sphere at W = 1:55, p) = 1 atmos- 
phere. The location of the fringe formation shows that for a supersonic 
projectile the disturbance is concentrated almost entirely in the region between 
two shock waves and in the narrow turbulent wake. 
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Figure 9. Interferogram of 20° included angle, calibre 0-225 cone-cylinder at 
M = 2:25, po = 1 atmosphere. 
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Figure 10. Interferogram of } in. sphere at MW = 1-35, py = 1 atmosphere. 
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D. H. Steininger and F. D. Bennett, Interferometric test for N-wave flow, Plate 4. 
: HH 
ee Figure 11. Interferogram of |} in. sphere at 17 — 1-29, p») = 1 atmosphere. 
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Yo =Xy—ar; (3) A from A = (x,—xy)/2r!4. We obtain an average 4 
and yy). Through a plot of the points (x,7) a straight line is drawn, anc 
the slope of this line is compared with « to make sure that any difference 
between them is less than the standard deviation of the «’s found by taking 
numerous pairs of points from the set which determines the line. 
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Figure 12. Shock positions, } in. sphere. 
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Figure 13. Shock positions, cone-cylinder. 
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Now the values yg, 4, and « are used in the shock equations 
= ort Ars, Xy = art+yo+ Ari, 
For each value of r for which measurements were obtained, x, and x, are 
calculated. We compare calculated x,, x, with measured xy, xg for the 


4 in. sphere in table 4. Figures 12 and 13 show the measured and calculated 
shock positions for the sphere and the cone-cylinder respectively. 


r Measured x | Measured x, 
2-230 0-7831 0-364 4-721 4-332 
3-676 1-719 1-416 6-222 5-912 
4-916 2-610 2°377 7-451 7-211 
6-276 3-666 3-463 8-789 8-601 
4-829 4-034 10-203 10-044 
9-626 6-385 6-223 12-041 11-941 

11-062 7-527 7-430 13-433 13-350 
12-552 8-789 8-691 14-793 14-803 
14-053 10-029 9-971 16-239 16:257 
15-532 P-225 17-686 17-683 
16-925 12-411 12-436 18-991 19-022 
18-415 13-694 13-723 20-460 20-449 
19-948 15-010 15-051 21-928 21-913 
21-384 16-316 16-298 23-233 23-280 
22-874 17-621 17-596 24-626 24-696 
24-723 19-154 19-218 26°355 26-458 
26-214 20-481 20-514 27-802 27-860 
28-030 22-167 22-106 29-629 29-576 
29-520 23-407 23-414 30-999 30-982 
31-141 24-854 24-839 32-555 32-509 
32-707 26-246 26-218 34-034 33-982 
34-350 27-562 27-666 35-503 35-526 
35-851 29-020 28-991 36-938 36-935 
37-014 30-010 30-018 38-015 38-026 
38-135 30-934 31-009 39-146 39-077 
39-821 $2°522 32-500 40-767 40-656 


Table 4. Comparison of N-wave shocks with actual shock positions for | in. sphere. 
All units are in projectile diameters. 


4.3. Fringe-shift measurement, +; in. sphere 

The interferogram (figure 11, plate +) is oriented on the ‘Telereader 
so that the x-axis of the coordinate system is along the projectile axis and 
the shock positions are measured according to the procedure of §4.2. The 
slope of the straight line through the points (x»,r) is used to obtain «. 
Actually, the drag on the ; in. sphere is so great that the effective Mach 
number along the length of the shock waves is not constant, and three 
values of x are determined: one for the upper third of the shock wave 
region (this value, « = 0-810, is given in table 3), one for the middle third, 
x = 0-800, and one for the lower third where x = 0-786. 
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‘lo interpret the interferogram the fringe shift 5 is measured along the 
path traced by a fringe between the shock waves. If r is the radial distance 
to the disturbed fringe at a certain point (x,7) and r, the distance to the 
same fringe in the undisturbed region in front of the shock, then 6 = (r —r,,)/A, 
where \ is the distance between the undisturbed fringes. ‘The coordinates 
(v, 7) along the fringe are measured for at least six values of x ranging from 
the front to the rear shock. For each value of x, six separate measurements 
are taken of the radial distance and averaged for r. Then 6 is calculated at 
cach point and multiplied by the scale factor. 


fringe trace 


Figure 14. Correction of fringe-shift measurement. 


‘lhe theory with which we wish to compare these measurements predicts 
tringe shift along a horizontal trace, not along a fringe; so we correct our 
measurements appropriately. It was shown in §2.5 that 6G is constant 
along the lines 

x = art Arl4(2é— 1). 


‘The slope of these lines deviates from « by an amount that ranges from 
Ar at the shocks to zero midway between the shocks where € = 0:5. 
For large r, Ar-** is small compared with «, and to a good approximation 
5G is constant along lines of slope «. We therefore move the measured 6G 
at some point on the fringe, say A(x,7r) in figure 14, along a line of slope « 
to a new point, say B(x’,r,,), on the horizontal trace. By this construction, 
v= x—2(r—r,). The value of € corresponding to the 5G measured at 4 
is then € = (x’—xy)/(xys—xXy). The change in r'® of the scale factor due 
to the move from point 4 to point B is less than 0-1°% for r > 70 diameters, 
and is neglected. Six separate measurements were taken of the distance 
between the shock waves (x, — xy), and their average was used in the calculation 
of €&. The value of 4 used in the scale factor is an average value determined 
from the shock wave positions in the manner described in § 4.2. 

The results of these fringe-shift measurements for several traces through 
the shock waves are given in figures 15 to 19, where 4G is plotted against ¢ 
and compared with the curve h(é) = €'7(1—3€). 
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Figure 15. Fringe shift between the shocks. 
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Figure 16. Fringe shift between the shocks. 
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Figure 18. Fringe shift between the shocks. 
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Figure 20. Fringe shift between the shocks. 
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4.4. Fringe-shift measurements, cone-cvlinder 
The interferogram of figure 9 (plate 2) was measured in the same manner 
as that of the sphere. The results are shown in figures 20 and 21. 
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Figure 21. Fringe shift between the shocks. 


5. DiIscUsSION OF RESULTS AND CONCLUSIONS 
5.1. Shock positions 

‘Table + compares the actual shock positions for the } in. sphere with 
the shock positions predicted for the \-wave by (2.2). Figure 12 presents 
this comparison pictorially. ‘lhe calculated points plot a smooth curve, 
but the measured points have a randomness due to errors in the measure- 
ments. At small 7, the measured points of the shocks fall well outside the 
N-wave curves of (2.2). For r > 15 diameters, the measured points fall 
randomly on either side of the N-wave curves and within ;, of a diameter 
from them, which is estimated to be within the probable error of the 
measurements. 

We assume the validity of (2.2) in order to determine the axis of the 
projectile and to measure the constants A and yy. We now conclude that 
this assumption is correct, and that the actual shock waves follow the curves 
given by (2.2) for all r > 15 diameters. 

Whitham (1952) has interpreted the experimental results of DuMond 
et al. (1946) as showing that the shock waves follow the curves of (2.2) for 
r > 1000 diameters. ‘The results given here indicate that good agreement 
is obtained closer to the projectile axis by a factor of 60. 
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A similar comparison of actual and .V-wave shock positions for the 
slender cone-cylinder is shown in figure 13. For this projectile the curves 
coincide for r > 7 diameters. 


5.2. Fringe shift in the N-wave 
Figures 15 to 19 show the results of the fringe-shift measurements along 
several traces of constant r, and compare these results with the fringe shift 
predicted by the theory as represented by the curve A(€) = €'7(1— 36). 
For small r, the measured data lie well above A(€) near the front shock and 
below A(€) near the rear shock. We attribute this to the fact that in the 
near regions where the ideal .V-wave has not yet formed, the pressure 
profile between the shocks lies above a linear profile near the front shock 
and below near “he rear shock (see figure 3). As 7 increases, the measured 
fringe-shift curve moves closer to A(é). At r = 74 diameters, the measured 
data fall very close to the exact theoretical curve of (2.8) (as represented by 
the dashed line of figure 15). The circles plotted in figure 15 are the average 
values of the measurements. ‘The lengths of vertical arrows through these 
circles represent the composite estimated errors in the fringe-shift measure- 
ments and in the determination of the scale factor. ‘The horizontal arrows 
show the estimated error in €. A detailed error estimate is given elsewhere 
(Steininger 1956). 

We notice that for r = 74 and r = 67, in the region where fringe shift 
all along the trace is rather small, the data points just behind the front shock 
fall below the theoretical curve. In the near regions of both cone-cylinder 
and sphere within 10 to 15 diameters of the axis, the data points do not fall 
below the theory. However, in determining the value of € for each fringe- 
shift measurement (see $4.3), we moved the value of fringe shift down a 
line of slope x to the horizontal trace, whereas to be completely accurate 
we should have moved down a line of steeper slope. ‘his means that a 
more accurate plotting of the data would move the points near the front 
shock wave to the right and hence below the theoretical curve. ‘This 
correction would increase with increasing fringe shift since it is proportional 
tor—r,, and with decreasing r since it is proportional to dr 4. An analysis 
of the errors involved in the measurements does not indicate the possibility 
of these low points. ‘Therefore we must look for a physical basis for the 
discrepancy. A possible explanation is that refraction etfects immediately 
behind the shock front displace the fringes slightly downward and result 
in a smaller fringe shift. Bennett e¢ a/. (1952) give evidence that fringe- 
shift data near the shock usually lie below expected values. More recent 
unpublished investigations indicate that refraction could cause measured 
values to lie either below or above expected values depending upon which 
side of the focal plane the projectile lies. ‘The position of our projectile 
with respect to the focal plane is not known, so we cannot come to a definite 
conclusion about the effects of refraction in our particular case. 

We conclude that (1) equation (2.8) does predict the fringe shift in the 
region between the shocks where the pressure profile is linear, (2) the 
approximate function h(é) = €'2(1— 3) can be used as a test for the 
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existence of .V-wave flow, and (3) for the sphere, .V-wave flow exists for r 
greater than about 70 diameters. We note in passing that the shock waves 
assume their asymptotic shapes much closer to the projectile axis than the 
pressure profile does. 

The results of fringe-shift measurements for the cone-cylinder are 
presented in figures 20 and 21 for qualitative comparison with the results 
for the sphere. We see that as r increases, the peak value of the measured 
fringe shift moves away from the theoretical curve, giving evidence that the 
actual fringe shifts, and consequently the pressures in this region, are 
dissipating less rapidly than an \-wave pressure distribution which the 
theory represents. This is to be expected. Whitham (1952) shows that 
the pressure curve WA(71), if it applied at small r, would become infinite 
at r = 0, and that the shock strength at large r, Wh(44), decreases at a rate 
proportional to r*4. It seems reasonable that, in the case of the cone- 
cylinder where the conical region of compressive flow extends far up the 
front shock wave, the actual pressure, although having a finite value at 
r = 0), dissipates much less rapidly at small values of r than the N-wave 
pressure would. 

In contrast, we notice that for the sphere where the flow is compressed 
at the front and then expands continuously around the body, the measured 
fringe shift converges upon the theoretical curve even in the relatively 
near regions. 


APPENDIX I. INTERSECTIONS OF THE FUNCTIONS f AND Rg 
From § 2.4, 


f=(e-1)" (1.1) 

and g = log[e + (e?—1)!?]. (1.2) 
The functions f and Rg intersect when 

(e?—1)'? = Rlog[e +(e? — 1)!?] (1.3) 

or exp[(e? — 1)!?/R] = «+ (1.4) 

Let 0 = (e?—1)!?/R. Then 

= €+(e2—1)!2, (1.5) 

and e~" = [e+ = e—(e2- 1)!2, (1.6) 

Now = = 112, (1.7) 

so sinhé = R@. (1.8) 


Figure 22 shows the relationship of the functions sinhé@ and R@. For the 
function sinh 


7p (sinh 6)>1 for@>O. 
For the function Ré, 
d(R?) 


| 
q 
7 
| 
| 
$ 


Interferometric test for N-wave flow 


In these two cases, the only intersection of R@ with sinh@ occurs at 
9=O(e=1). When R>1, 


and R@ and sinh @ intersect twice, once at 6 = 0 and once for 0 > 0 (e > 1). 
A 


Figure 22. Relationship of the functions sinh @ and R@. 


APPENDIX II. VELOCITY CALCULATION FOR SPHERES 
‘The equation of motion for any non-yawing projectile in free flight can 


be written 
dV 


dz 
where V is the velocity of the projectile, and A, the ballistic coefficient ; 
also Fy = po d?/m, where pg is the air density, d the diameter of the projectile, 
and m the mass of the projectile. From the work of Charters & Thomas 
(1945), for spheres, 
Ky = (0:3812 + 0-0006) — (0-0140 + 0-0008)( — 2-75), 
or in general terms Kp = a(1—5dV), (11.2) 
where a = 0-4197 and 6 is a constant which depends on the sound velocity. 
Equation (II.1) now becomes 

= —F,aV(1-dV), (11.3) 

which can be integrated to give 
= = V, exp( — F, az). (11.4) 


(1.1) 
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We replace |’ by dz/dt, integrate again and find that 
ex p( Fy az) 
Fya 
where C is the second constant of integration. 

For the sphere shown in figure 10 (plate 3), the time ¢, required to 
traverse a distance s, between the two velocity stations of the range was 
measured. We evaluate the constant C in (11.5) from the initial condition 
that when ¢ = 0, s = 0, and the constant V’, from the final condition that 

Now the distance 3, from the first velocity station to the interferometer 
is measured, and from (II.4) the value of V at the interferometer is determined 
thus 


+bV,z = V,t+C, (11.3) 


A 
= V, exp( — Fy (11.6) 
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SUMMARY 


‘his paper is concerned with the problem of obtaining higher 
approximations to the flow past a sphere and a circular cylinder 
than those represented by the well-known solutions of Stokes 
and Oseen. Since the perturbation theory arising from the 
consideration of small non-zero Reynolds numbers is a singular 
one, the problem is largely that of devising suitable techniques for 
taking this singularity into account when expanding the solution 
for small Reynolds numbers. 

‘lhe technique adopted is as follows. Separate, locally valid 
(in general), expansions of the stream function are developed for 
the regions close to, and far from, the obstacle. Reasons are 
presented for believing that these ‘ Stokes’ and ‘ Oseen’ expansions 
are, respectively, of the forms 


9) and > (Rr, 4) 

where (r,@) are spherical or cylindrical polar coordinates made 
dimensionless with the radius of the obstacle, R is the Reynolds 
number, and f,,.,/f, and F,,.,/F,, vanish with R. Substitution 
of these expansions in the Navier-Stokes equation then yields a 
set of differential equations for the coefficients ys, and ‘V’,, but 
only one set of physical boundary conditions is applicable to each 
expansion (the no-slip conditions for the Stokes expansion, and 
the uniform-stream condition for the Oseen expansion) so that 
unique solutions cannot be derived immediately. However, the 
fact that the two expansions are (in principle) both derived from 
the same exact solution leads to a ‘matching’ procedure which 
yields further boundary conditions for each expansion. It is thus 
possible to determine alternately successive terms in each 
expansion. 

The leading terms of the expansions are shown to be closely 
related to the original solutions of Stokes and Oseen, and detailed 
results for some further terms are obtained. 


1. INTRODUCTION 

‘The problem of determining the steady flow past fixed bodies in a slow 
uniform stream of viscous incompressible fluid is an old one. It was first 
considered by Stokes (1851), and has been discussed subsequently by 
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many authors. With very few exceptions, however, these authors have 
been almost entirely concerned with finding the flow past various shapes 
of body in the limit of zero Reynolds number. Yet many of the effects 
that arise when the Reynolds number is not negligibly small are also of 
considerable physical and mathematical interest. A number of phenomena 
in lubrication and the motion of small particles, for instance, depend 
critically on second-order effects arising from the inertia of the fluid. For 
practical purposes, a first approximation to these second-order effects 
would doubtless be quite adequate, but the purely mathematical difficulties 
encountered in their calculation raise fundamental questions concerning 
the general nature of expansions for flow fields at small Reynolds numbers. 
It is with such questions that the present paper is concerned. 

The problem originally considered by Stokes (1851) was that of flow 
past a sphere, for which he obtained a solution by neglecting completely 
the inertia of the fluid. Later, Whitehead (1889) attempted to improve 
upon this solution by obtaining higher approximations to the flow when 
the Reynolds number is not negligibly small. The method proposed by 
Whitehead was the obvious one of using a lower-order approximation to 
calculate the inertia terms in the equation of motion, thus developing an 
iterative procedure. Since the boundary conditions at each stage of the 
iteration are independent of the Reynolds number, this procedure is clearly 
equivalent to assuming an expansion of the flow in powers of the Reynolds 
number. When this assumption is valid, and there are many slow motion 
problems for which it is*, there is little more to be said. But, as is now 
well known, the assumption is never valid in problems of uniform streaming. 
‘The particular difficulty encountered by Whitehead was that the second 
approximation to the velocity of flow past a sphere remains finite at infinity 
in a way which is incompatible with the uniform-stream condition. And 
higher approximations to the velocity distribution actually diverge at 
infinity. ‘The assumption of an expansion in powers of the Reynolds number 
thus leads to a situation in which it is impossible to satisfy the boundary 
conditions of the problem in all terms except the leading one. This 
mathematical phenomenon appears to be common to all problems of 
uniform streaming past bodies of finite length-scale, and is sometimes 
referred to as ‘ Whitehead’s paradox’. 

The paradox has, of course, long since been resolved. Both its physical 
origin and a mathematical device for overcoming the associated difficulties 
were pointed out by Oseen (1910). Since any body moving steadily through 
a viscous fluid must experience some resistance, consideration of the 
momentum flux across a large surface surrounding the body shows that 
the magnitude of the disturbance to the uniform velocity of the stream 


* Apart from some finer qualifications, the necessary and sufficient condition for 
the validity of this assumption in problems of steady flow is that the velocity should 
fall to zero at a great distance from the boundaries generating the flow not less rapidly 
than the reciprocal of 'the distance. ‘The assumption is thus valid for the flow 
generated by a slowly rotating sphere. 
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cannot everywhere fall to zero more rapidly than the inverse square of the 
distance from the body. But, for so large a disturbance as this, the 
acceleration of the fluid at a great distance cannot be negligible by comparison 
with the local viscous force. For the former, being almost entirely due to 
the convective effect of the stream, is a constant multiple of the first 
derivative of the velocity, whereas the latter is a multiple of the second 
derivative of the velocity. ‘Thus the viscous force can be dominant only 
if the decay of the disturbance is exponentially rapid. Stokes’s theory is 
therefore not self-consistent at a great distance from the body, and it is not 
surprising that the procedure adopted by Whitehead should lead to further 
inconsistencies. In more mathematical terms, the perturbation represented 
by asmall non-zero Reynolds number has a singularity at infinity (in space), 
and Stokes’s solution does not provide a uniformly valid approximation to 
all the required properties of the flow. 

It is a straightforward matter to show that Stokes’s solution does not 
break down until the region in which the flow is nearly a uniform stream 
has been reached, so that the solution does provide a uniformly valid 
approximation to the total velocity distribution (and consequently a valid 
approximation to many bulk properties of the flow, such as the resistance). 
It is only the derivatives of the velocity at a great distance that are seriously 
in error. But, of course, this error is crucial in the problem of obtaining a 
second approximation to the flow, since the neglected inertia terms in the 
equations of motion involve velocity gradients. Nor can Whitehead’s 
procedure be used to derive a locally valid second approximation to the 
flow in the region not far from the sphere. It is true that this procedure 
produces a correct differential equation for the second approximation, but 
the spatial restrictions placed upon its validity prohibit the use of the outer 
boundary conditions, so that a unique solution cannot be derived. In fact, 
as in all singular perturbation problems, a uniformly valid approximation 
to the neglected terms in the governing equation is a necessary prerequisite 
for the determination of a second approximation to the solution anywhere 
in the field. 

Fortunately, as was shown by Oseen (1910), the determination of a 
uniformly valid first approximation to the velocity and all its derivatives 
is itself a linear problem which may be solved analytically. ‘The circum- 
stance, already mentioned, that the inertia terms are important only in the 
region where uniform-stream conditions have been almost attained permits 
a linear approximation to be made which yields the well-known Oseen 
equation. he very interesting description of the asymptotic flow field 
given by the solutions of this equation now occupies an important place 
in the theory of viscous motion. In the present paper, however, we are 
more interested in the use of these solutions as tools in the problem of 
finding higher approximations. 

We may also note that the relevant solution of Oseen’s equation provides 
a uniformly valid approximation to the velocity and all its derivatives in 
the two-dimensional flow past an infinite cylinder of finite cross-sectional 
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length-scale, though Stokes’s approximation yields no solution in this case. 
‘he first such solution to be obtained was that for flow past a circular 
evlinder (Lamb 1911), and a more quantitative summary of these results 
and those for flow past a sphere is given in § 2. 

In both the two- and three-dimensional cases, therefore, the solutions 
ot Oseen’s equation provide an adequate starting point for the determination 
of higher approximations to the flow. However, apart from the recent 
paper by Lagerstrom & Cole (1955), no work on this problem seems to 
have been published. It is true that a great deal of effort has been expended 
in finding higher approximations on the assumption that Oseen’s equation, 
rather than the Navier-Stokes equation, is the exact governing equation ; 
but this problem is not strictly relevant. Oseen’s equation does contain 
the Reynolds number as a free parameter (an obviously necessary conse- 
quence of the uniform validity of Oseen’s approximation), but the idea of 
solving the equation to a higher order of approximation in the small Reynolds 
number than that involved in its derivation is of limited value. ‘The 
justification usually given is that Oseen’s equation and the Navier-Stokes 
equation are qualitatively similar, so that solutions of the former might be 
expected to yield qualitative information about solutions of the latter for 
all Reynolds numbers. And on these grounds, a few basic solutions, such 
as Goldstein’s (1929) solution for the sphere, are surely worth while. But 
the problem as a whole has probably received far more attention than it 
deserves. 

In principle, the problem of obtaining higher approximations to the 
real flow is not appreciably more difficult than that mentioned above. 
For there seems little reason to doubt that Whitehead’s iterative method, 
using Oseen’s equation rather than Stokes’s equation would yield 
an expansion, each successive term of which would represent a 
uniformly valid higher approximation to the flow. In each step of the 
iteration a lower-order approximation would be used to calculate those 
particular inertia terms that are neglected in Oseen’s equation and the 
resulting inhomogeneous form of Oseen’s equation would be solved, to 
the relevant degree of accuracy, for the boundary conditions appropriate 
to the problem. ‘The expansion generated in this way would seem to be 
the most economic expansion possible, in the sense that the partial sums 
of any order contain all the legitimate, and no redundant, information about 
the whole flow field. 

Nevertheless, it is the primary object of the present paper to describe 
in detail an alternative procedure which in many ways is more satisfactory. 
This alternative procedure involves simultaneous consideration of locally 
valid (in general) expansions close to, and far from, the singularity of the 
perturbation. ‘These expansions may be called ‘Oseen’ and ‘Stokes’ 
expansions, respectively, since their leading terms are closely related to the 
original approximations of these authors. The Stokes expansion is a 
straightforward expansion of the kind envizaged by Whitehead. It is an 
expansion in the Reynolds number for fixed values of the space coordinates 
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{made dimensionless with the finite length-scale of the body). For the 
Oseen expansion, on the other hand, the coordinate system is first strained, 
by a factor depending upon the Reynolds number, in such a way that the 
length-scale of variations in the asymptotic flow pattern at a great distance 
from the body is finite in terms of the new coordinates. In this new 
coordinate system, the length-scale of the body is very small, and the 
singularity of the perturbation is removed to the origin of coordinates 
(inside the body). ‘The Oseen expansion is then an expansion in the 
Reynolds number of fixed values of the new coordinates. The connection 
between this expansion and Oseen’s work is evident when one remembers 
that both the inertia and viscous forces at a great distance depend J/inearly 
on the disturbance to the stream. For the fact that these forces are of a 
comparable order of magnitude is then made to appear ‘natural’ by the 
choice of length-scale referred to above. 

On grounds of expediency alone, the use of Stokes and Oseen expansions 
is preferable to the use of an expansion which is generated by uniformly 
valid successive approximations. In the first place, their mathematical 
structure is a great deal simpler. In fact, they are usually power series, or 
simple extensions of power series, in the Reynolds number; whereas a 
uniformly valid approximation necessarily depends upon the Reynolds 
number in a complicated manner since it involves functions of both the 
‘strained’ and ‘unstrained’ coordinate systems. Moreover, uniformly 
valid approximations per se are not usually of much physical interest (apart, 
perhaps, from a first approximation, which gives an overall description 
of the singularity of the perturbation). In the present problem, for instance, 
it is the Stokes expansion that gives virtually all the physically interesting 
information. Questions of uniformity arise only in connection with the 
proper derivation of this expansion, and there is therefore some point in 
attempting to cast such questions——and techniques for answering them— 
in terms of the expansion itself. 

However, the attitude adopted by Lagerstrom & Cole (1955), who also 
discuss in general terms the use of Stokes and Oseen expansions in the 
problem of the present paper, is rather more fundamental. They point 
out that, when dealing with asymptotic expansions for small Reynolds 
numbers, it is wise to restrict attention to those expansions that can (in 
principle) be derived from the exact solution by the application of formal 
limit processes which may be defined @ prior’. For it is then a relatively 
straightforward matter to discuss such questions as the error involved in 
any particular partial sum, or the domain of uniform validity of the 
expansion. The Stokes and Oseen expansions are of this type since they 
may be derived by the limiting processes described above. ‘Thus, 
if R is the Reynolds number and w is the (dimensionless) position vector, 
the limiting process R->+0 for either fixed x or fixed Ry defines, 
respectively, the Stokes or Oseen limit. The application of these limits 
to the difference between the exact solution and the wth partial sum then 
yields the (7+ 1)th term of the relevant expansion. 
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An expansion in terms of uniformly valid successive approximations, 
on the other hand, cannot normally be derived from the exact solution in 
this way. It can be derived from the exact solution by inspection, when 
the detailed form of the solution is known; or it can be defined, as has 
been done above, in terms of the iterative process that is to be applied to 
the governing equation. However, both these concepts are very cumbersome, 
and since uniformly valid approximations may always be constructed from 
a simultaneous knowledge of the Stokes and Oseen expansions, it is more 
satisfactory to proceed in this latter way rather than vice versa. 

The problem thus reduces to that of determining the proper boundary 
conditions for the individual terms of the Stokes and Oseen expansions. 
At the regular end of the range of validity of these expansions (the body 
for the Stokes expansion, and infinity for the Oseen expansion), the 
boundary conditions are the physically obvious ones. At the singular 
end, however, the physical boundary conditions are irrelevant and it is 
necessary to use the fact that the Stokes and Oseen expansions are different 
forms of the same function. ‘This leads to a matching of the expansions 
which is of such a kind that it becomes possible to derive alternately successive 
terms in each expansion. ‘These matching conditions have already been 
described in general terms by Lagerstrom & Cole (1955) and a detailed 
account of the procedure is now given in §3 and § 4. 

For the sake of simplicity, the paper deals only with flow past a sphere 
and a circular cylinder (treated respectively in $3 and $4), since these special 
cases appear to illustrate most* of the pertinent ideas. 


2. 'THE APPROXIMATIONS OF STOKES AND OSEEN 


2.1. Flow past a sphere 

In Stokes’s original treatment (1851) of slow streaming past a sphere, 
the inertia of the fluid is neglected completely, so that the Navier-Stokes 
equation reduces to 

0 = —gradp+vV7u, (2.1) 

where p is the kinematic pressure, v is the kinematic viscosity, and u is the 
velocity vector. It is then a straightforward matter to obtain, in terms of 
a Stokes stream function, the integral of this equation that satisfies the 
no-slip condition at the sphere and the uniform-stream condition at infinity. 
If U is the velocity of the undisturbed stream, and a is the radius of the 
sphere, this integral for the Stokes stream function Ua is 


1 1 
3r+ -)sin @, (2.2) 
where the origin of polar coordinates (ar, 6, ¢) is taken at the centre of 


the sphere, the line 6 = 0 being in the direction of the undisturbed stream. 


*'Though certiinly not all. For a recent account of the important effects of 
asymmetry in two-dimensional flow past a cylinder, see Imai (1951). The corre- 
sponding three-dimensional problem appears not to have been considered. 
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The individual terms in (2.2) are, successively, a uniform stream, a so-called 
‘stokeslet’, and a doublet. It is, of course, only the stokeslet which 
contributes to the vorticity of the flow. 

The values of r at which Stokes’s approximation breaks down may 
now be found by using (2.2) to calculate the neglected terms. When + is 
large, the dominant inertia terms are the convective accelerations arising 
from the interaction between the uniform stream and the stokeslet, and 
are of order U?/ar?. The viscous term, on the other hand, is of order 
vU/ar* in this part of the flow. The Stokesian theory therefore breaks 
down when 

Rr = O(1), (2.3) 


where R is the small Reynolds number Ua,v. It is important to note in 
this criticism that the dominant inertia terms arise from the rotational 
velocity field due to the stokeslet, and cannot be represented as the gradient 
of a scalar. If the latter had been possible, the error arising from these 
terms could have been absorbed into the solution for the pressure field 
without affecting the solution (2.2) which is the general solution for a 
conservative distribution of viscous forces. 

According to (2.3), the solution (2.2) approaches arbitrarily close (for 
sufficiently small R) to the conditions of a uniform stream before the theory 
breaks down. Stokes’s solution is therefore actually a uniform approxi- 
mation to the total velocity distribution*. However, it is clearly not a 
uniform approximation to the disturbance of the uniform stream, or, 
equivalently, to the distribution of velocity gradients. As noted in §1, 
it is for this reason that the solution cannot be used to obtain a second 
approximation in the manner attempted by Whitehead (1889). 

The idea behind Oseen’s (1910) technique for obtaining a uniform 
approximation to the disturbance of the stream is to take inertia forces 
into account in the region where they are comparable with viscous forces, 
but neglect them in the Stokesian region of the flow. Thus, since the flow 
is nearly a uniform stream in the former region, the appropriate equation is 


U.gradu = —gradp+vV7u, (2.4) 


where the vector U represents the uniform stream. ‘lhe left-hand side 
of (2.4) is, of course, negligible throughout the region in which Stokes’s 
approximation is valid. It may be noted in passing that the equation (2.4) 
is formally the same as the equation which would be obtained if the velocity 
distribution were written in the form U +u and the Navier-Stokes equation 
were linearized in the disturbance velocity u. However, this interpretation 
is conceptually wrong and can lead to erroneous or misleading conclusions 
such as Lamb’s statement (1932, p. 610) that Oseen’s theory is less accurate 


* In this paper, we consider u* to be a uniform approximation to u if, as R -- 0, 
u—u* = 0 (u) for all values of x. This condition may clearly be violated in a trivial 
sense near the regular zeros of u. However, in any particular problem, we are 
interested in the condition only near the singularities of the perturbation. 
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than Stokes’s in the neighbourhood of the sphere (where the boundary 
condition u = —U would make nonsense of such a linearization). The 
left-hand side of (2.4) is not intended to be a uniform approximation to the 
inertia terms, and the difference between Oseen’s and Stokes’s theory in 
the neighbourhood of the sphere is of a small order which neither theory 
is entitled to discuss. 

The derivation of the appropriate exact solution of Oseen’s equation (2.+) 
is a matter of some difficulty (see Goldstein 1929). Fortunately, however, 
there is no justification (in the present investigation) for finding a solution 
which satisfies the boundary conditions to a higher order of approximation 
than that involved in the governing equation, and it is a relatively simple 
matter to show that the solution given by Oseen himself is an adequate 
approximation. In terms of the dimensionless stream function yf, this 
solution is 

‘The function (2.5) certainly satisfies the differential equation (2.4) and the 
relevant boundary condition at infinity. Moreover, when r is of order 
unity, the function becomes 


1 
4(2°-3r+ )sin*0 + O(R), 


which agrees with Stokes’s solution, and consequently satisfies the boundary 
condition on the sphere, to an adequate approximation. 

That (2.5) provides a uniform approximation to the disturbance of the 
stream follows immediately from the manner in which it has been obtained, 
and it is interesting to observe the nature of the non-uniformity associated 
with Stokes’s solution. According to (2.5), the symmetric flow due to the 
stokeslet changes, at large values of Rr, into a simple source described bv 


the stream function 3 
brn — 7p + cos 8), 


whose mass fl»x is supplied asymmetrically by an inflow along a narrow 
wake defined, in an order-of-magnitude sense, by 


Rr(1—cos @) = O(1). 


‘The function that describes this transition is a function of Rr and 6, and 
this strained coordinate system will appear as the natural one for the entire 
Oseen expansion (see $3.2). 

Finally, it should be noted that Oseen’s equation (2.4) is only a first 
approximation to the governing equation and cannot be used as such to 
derive second approximations to any property of the flow. Thus Oseen’s 
derivation (1913) of a second approximation to the drag coefficient of a 
sphere, namely 


Cy = Tap = sk), (2.5) 
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where D is the kinematic drag, would seem to be invalid, though the result 
is in fact correct (see $3.4). Strictly, Oseen’s method gives only the leading 
term of (2.6) and is scarcely to be counted as superior to Stokes’s method 
for the purpose of obtaining the drag. 


2.2. Flow past a circular cylinder 

The two-dimensional case exhibits profound differences. For slow 
uniform streaming motion past a circular cylinder, there is no solution 
to Stokes’s equation (2.1) that remains finite as r becomes indefinitely large 
and that satisfies the no-slip condition on the cylinder. Hence, unlike the 
three-dimensional case, there is no solution to Stokes’s equation that 
provides a uniform approximation to the total velocity distribution; it 
is therefore not immediately clear that Oseen’s equation (2.4) will provide 
a uniform approximation to even the total velocity distribution. 

Let ar and 6 refer in this section to the two-dimensional radial and 
angular polar coordinates respectively, the line 6 = 0 being the positive 
direction of the stream, and let Ua’ be the Lagrange stream function, 
where a is the radius of the cylinder and U the uniform streaming velocity. 
If, now, we seek that solution of Stokes’s equation which satisfies the 
no-slip condition on the cylinder r = 1, which contains a term sin@ (in 
view of the uniform stream condition at infinity), and which diverges least 
rapidly as r becomes indefinitely large, we obtain 


f= c(; —r+2rlog r) sin 0. (2.7) 


For very large r, the solution (2.7) is dominated by the rotational term 
rlogr, and, at first sight, appears wholly unsatisfactory because of this 
logarithmic divergence. However, by substituting (2.7) into the full 
equation of motion, we find that the inertial forces neglected are of order 
C*(log r)/r?, and the viscous forces of order C/Rr®. These will be comparable 
when 


CRr logr = O(1), (2.8) 


and we should not expect the Stokes solution (2.7) to be valid beyond 
a value of r given by (2.8). In other words, while (2.7) may be an adequate 
representation of the flow relatively close to the cylinder, it cannot represent 
a uniform approximation to the total velocity distribution. 

It is possible, however, to write the approximation (2.7) in such a way 
that this non-uniformity appears less severe than might at first sight be 
supposed. ‘Thus we may write (2.7) in the form 


p= c( — 2rlog{f(R)} + 2rlog{rf(R)} —r+ ~)sin 6, (2.9) 


when f(R) is an arbitrary function of R. For f(R) < 1, and 7f(R) of order 
unity, the dominant term in (2.9) will be —2Clog{f(R)\r sin@. If this 
is to represent the external flow, namely a uniform stream rsin@, then 
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we can put C = —1/2log{f(R)}. Further, by substituting this value of ( 
and r = 1/f(R) into (2.8) we get 
R/f(R) = O(1). (2.10) 
Thus for f(R) = R, the Stokes solution (2.9) leads to a uniform stream 
of order unity in that very region where Stokes’s equation ceases to be 
valid. It is this apparently fortuitous behaviour that suggests that the 
external uniform stream condition has been reached before the Stokes 
approximation breaks down. ‘The logarithmic term in (2.7) apparently 
plays an important role in making the solution (2.7) as nearly uniform as 
the approximation allows, and it is interesting to note that a uniform 
approximation is obtained by making only a slight change in the solution 
(2.7). Thus the form 
1 1 r ‘ 
— | - sin @ 
is such that it tends to r sin @ as r tends to infinity, and also only differs from 
the Stokes solution (2.7) by terms of order R when r is of order unity. 
The form (2.11) may be regarded as an analogue of (2.2) in the three- 
dimensional case, and for similar reasons, we are led to expect that Oseen’s 
equation (2.4) will provide a uniform approximation to the disturbance 
stream function. 
In terms of the stream function, Oseen’s equation becomes 
Cx 


where x =rcos@. Faxén (1927) and later ‘Tomotika & Aoi (1950) have 
solved this equation exactly to satisfy the no-slip condition on the cylinder 
and the uniform-stream condition at infinity. However, as has been 
explained for the case of the sphere, there is no point in solving the linear 
equation (2.12) to a greater degree of approximation than that of the inertial 
terms neglected by substituting the Oseen equations for the Navier-Stokes 
equations, and so the simpler solution given by Lamb (1911) is as good an 
approximation as it is possible to obtain from the linearized equation. 
In fact, by writing Lamb’s solution in terms of #/, in the form 


1 rsinné 
ab (: 2 Op, (2.13) 
where By = 4-—y—log{R (2.14) 
and $, = 4K, 1, +2K (2.15) 


the K,, and J, being modified Bessel functions, we get a uniform approxi- 
mation to the disturbance stream function over the entire flow. Within 
the Stokes region (2.13) becomes, to the approximation involved, a Stokes 
solution, while far from the cylinder it represents a characteristic wake flow 
similar to that described for the sphere. But whereas the solution for the 
sphere is correct to terms of order R, the solution for the cylinder is correct 
only to terms of order (log R)-}. 
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3. GENERAL EXPANSIONS FOR FLOW PAST A SPHERE 

In this section we consider the problem of how to obtain higher 
approximations to the velocity distribution in flow past a sphere. The 
disadvantages of applying a uniformly valid iterative process to Oseen’s 
equation (2.4) have already been noted in $1. The terms of the resulting 
expansion, of which the leading term (2.5) is typical, involve the Reynolds 
number in far too heterogeneous a manner to make the significance of the 
expansion at all clear. Instead, therefore, we employ the more homogeneous 
Stokes and Oseen expansions are described below. 


3.1. The Stokes and Oseen expansions 
In the Stokes region of the flow, that is where r = O(1), we write the 
governing equation for the stream function ¢ in its customary form, namely 


1 2. 1 
= cosd, (3.2) 
cr 
3 
and assume an expansion of the form 
= fo + A(R), + (3.5) 
where + (R) 'f,(R) 0 as R 0. (3.6) 


The expansion (3.5) should be regarded as the expansion of the exact 
solution ¢(r, 4, R) for small values of R at a fived value of r. The assumption 
that the expansion takes the form indicated is thus equivalent to the very 
mild assumption that there is no singular dependence on R in the finite 
part of the field (such as arises, for instance, in connection with the formation 
of shear layers at Jarge values of R). We refer to (3.5) as the Stokes 
expansion. 

Since we require the magnitude of the velocity to be everywhere bounded, 


we may write 
1 3.7) 


without loss of generality. Formal allowance is thus made for the possibility 
that 4, = 0, though in fact this is not necessary, since it is evident that 
us, must be Stokes’s solution. 

The expansion (3.5) is required to satisfy the differential equation (3.1) 
and the no-slip condition on the sphere. Since the expansion is invalid 
at large values of r, the uniform-stream condition at infinity must be replaced 
by the requirement that the expansion should be perfectly matched to an 
expansion which 7s valid in the outer region. 


oe 
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The reason for the breakdown of the Stokesian theory is that inertia 
and viscous forces become comparable at large values of 7. This suggests 
that, for the outer expansion, we should find a transformation which removes 
the Reynolds number from the governing equation, thereby suiting the 
coordinate system to the fact that all terms in the equation are of a comparable 
order of magnitude. ‘There are, of course, many such transformations, but 
for the purpose in hand we assume that they are all essentially equivalent 
and select the simplest, which is an isotropic straining of the coordinate 
system and a simple scaling of the stream function. ‘Thus we introduce the 
variables 

p=f(Ry, = (3.8) 
and, when the governing equation is expressed in terms of these variables, 
the condition that the R should not appear is 


Rf(R) = @(R). (3.9) 
Further, since the (dimensionless) velocity must be of order unity in the 
region of validity of the expansion now sought, that is, where p and ‘are 
of order unity, we must have 


F(R) = g(R). (3.10) 
In this way we obtain the Oseen variables 


in terms of which the governing equation (3.1) becomes 
1 a(¥, D5 
O(p, 
where Dz and L, are the same operators as (3.3) and (3.4), but with r replaced 
by p. 
‘The expansion in the outer region, which we call the Oseen expansion, 
is now assumed to take the form 


where as (3.14) 


+ = (3.12) 


‘That the leading term should be independent of R is, of course, implicit in 
the choice of the Oseen variables (3.11). The structure of the remainder 
of the expansion then rests upon assumptions similar to those made for the 
Stokes expansion. In the present case, however, the absence of a singular 
dependence on p and pu (except at p = 0) as R ~ 0 is not intuitively obvious, 
and the plausibility of the expansion (3.13) depends on the demonstration 
that it is possible to satisfy the governing equations and boundary 
conditions with such an expansion. 

‘The expansion (3.13) is required to satisfy the differential equation (3.12) 
and the uniform-stream condition at infinity. On the other hand, the 
no-slip condition on the sphere, which in the new coordinate system has 
shrunk to a very small sphere of radius R (thus introducing the Reynolds 


2% | 
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number into the boundary conditions rather than the differential equation), 
is replaced by the condition that (3.13) should be matched to the Stokes 
expansion (3.5) at small values of p. 

The matching conditions follow from the fact that the Stokes and Oseen 
expansions must be related in some way, since they are both expansions of 
the same function for small values of R, even though one expansion does 
not in general determine the other uniquely. The nature of this relation, 
and hence the conditions themselves, may be found from the following 
rather intuitive argument. The common features of the Stokes expansions 
of the infinitely many functions that all have the same Oseen expansion 
ought to be possessed by the ‘ Stokes expansion’ of the Oseen expansion 
itself. Thus, if the Oseen expansion 

is formally expanded about R = 0) for fixed values of r (by expanding the 
‘t'.(p, 4) for small values of p, and re-arranging the terms), the resulting 


expansion, g,(R)¢,,(r, Lt), 

must be closely related to the Stokes expansion (3.5) of the actual flow. 
In fact, the only reason why these last two expansions might not be identical 
is that the Stokes coefficients 4) might involve terms (like e~’ = e *) 
which are important in the Stokes region but which are so small when r is 
large that they do not contribute to any term of the Oseen expansion. 
‘Thus, _ all such possibilities into account, it seems that we must have 
f,(R) = g,(R) and that u,,(7, ~) and ¢,(7, 4) must have the same asymptotic 
expadalian for large values of r. ‘The expansions of the Oseen coefficients 
for small values of p thus determine uniquely the expansions of the Stokes 
coefficients for large values of r (and vice versa). 

In practice, of course, it is necessary to solve for one term at a time, 
and the procedure for, say, a Stokes term is then as follows. ‘The general 
solution for this term that satisfies the no-slip condition on the sphere is 
first obtained from the relevant differential equation and is then expanded 
as an Oseen expansion. ‘The individual terms in this expansion are then 
compared with the corresponding terms that have previously been calculated 
in the Oseen expansion of the full solution. According to the conditions 
derived above, the former terms must occur explicitly in the expansions 
of the latter terms for small values of p. A comparison of coefficients then 
uniquely determines the arbitrary constants in the general solution for the 
Stokes term. ‘The analogous procedure for the calculation of an Oseen 
term is obvious. 


3.2. The leading terms of the expansions 
It is clear from the conventional approach described in $2.1 that the 
leading terms of the expansions (3.5) and (3.13) are, respectively, Stokes’s 


solution 1 


: 
ang 
regs 
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and a uniform stream 

‘Fo = 2p°(1 (3.16) 
‘hese expressions both satisfy their respective ditterential equations and 
boundary conditions. Moreover, the matching of these leading terms is 
particularly trivial. When (3.15) is written in terms of the Oseen variables 
(3.11), its contribution to ‘Vis seen to be 


BY 
+ —p), (3.17) 


so that only the uniform-stream term contributes to I’). ‘Vhe boundary 
conditions for ‘I, are therefore the same at p = 0 and p = x, and the exact 
solution of the non-linear equation for ‘’y is the uniform stream (3.16). 
This result is, of course, merely a re-statement of the fact that Stokes’s 
solution gives a uniform approximation to the total velocity distribution. 

Although the forms of the leading terms are already known, it is 
instructive to consider their derivation ab initio, if for no other reason than 
to prepare the ground for the rather more difficult two-dimensional case. 
In this connection it should be noted that it is sufficient to determine either 
of the solutions (3.15) and (3.16), since adequate boundary conditions are 
then available for the determination of the other. 

Now, there appears to be a simple physical argument for the solution 
(3.16). The Oseen variable p is, by definition, 


where r’ is the (dimensional) radial distance from the centre of the sphere ; 
and this is independent of a. For fixed values of U and rv, therefore, a 
fixed value of p corresponds to a fixed position in space. Hence, if the 
limiting process R-~0 is interpreted as the limiting process a— 0, it 
follows that the flow at the fixed point under consideration must ultimately 
be that of the undisturbed stream. It is true that the argument makes some 
assumptions about the magnitude of the disturbance caused by a body of 
zero size, implying, for instance, that the total force on the sphere vanishes 
with its radius, but such assumptions can be accepted with some assurance. 
The derivation of the leading term of the Stokes expansion, x, is then 
straightforward. ‘lhe equation for uy is 
Dt by = 0, (3.19) 
and it is not difficult to show that the general regular solution of this equation, 
that vanishes at « = +1 has a double zero at r = 1, is 


be > [A,, \(2n = 1)r” (2n 1)r" 2r-n+2) 
+B, —(2n + *? + (2n— (3.29) 


where 4, and B, are constants, and 


O,(u) = 


P,(u) dp, (3.21) 


: 
| 
10 
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in which P,,(u) is the Legendre polynomial of degree n. When (3.20) 
is expressed in terms of the Oseen variable p, its contribution to ‘Vis found 
to be 

n=1 

+ B,{2R-"+1p"*1 — (2n+ 1)R"p-"*? + JO, (u), (3.22) 

and the requirement that this contribution should not contain terms of 
greater order than unity then gives 


A, «0, 
3.23) 
B, =0 (n > 2), i 
Further, since O,(u) = —3(1—p), (3.24) 


the requirement that the term of order unity should represent the uniform 
stream (3.16) gives 
B, = =}. (3.25) 


Thus, Stokes’s solution (3.15) is recovered. 


3.3. The second term of the Oseen expansion 

Since the leading term of the Oseen expansion is a uniform stream, the 
equation for the second term ‘’, must, in effect, be Oseen’s equation (2.4). 
In terms of the stream function, this equation is 

(3.26) 
p Cu a p 

It may be noted that this equation now appears naturally as the equation 
for the perturbation of a uniform stream, which appears to be at variance 
with the view taken in §2.1. The reason is, of course, that we are no longer 
dealing with a technique for obtaining a uniform approximation to the 
solution. 

Equation (3.26) may be solved by the method used by Goldstein (1929). 
The transformation 


= (3.27) 
reduces the equation to 
(D3 — = 0, (3.28) 
and the general solution that vanishes at infinity and « = +1 (the latter 
condition follows immediately from (3.3) of Y’ vanishes at u = +1) is 
easily found to be 
D 
(3.29) 


n=1 


where K,,,(4p) is a modified Bessel function. Since the order of the 
Bessel function is half-integral, it may be expressed in the closed form, 


(n+s)! 


) (n—s)!p*° 


.4(4p) = (3.30) 


9) 
d 
is 
\ 
) 
n 
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The completion of the general integration is rather troublesome (see 
Goldstein 1929). However, this is not necessary since it is possible to 
apply further boundary conditions in the partially integrated form (3.29). 
When (3.29) is expressed in terms of the Stokes variable r, its contribution 
to D*us is seen to be 

4, (ARNO, (1). (3.31) 
n=1 
Then, since (3.30) shows that 
(4p)!K,, ~ /p" (3.32) 
for small values of p, the condition that (3.31) should not contain terms of 
greater order than unity yields the results 


FAR) = R, (3.33) 
4, =0 (n > 2). (3.34) 
The contribution to D?% is thus 
"(1 —p?). (3.35) 
From Stokes’s solution (3.15), on the other hand, we obtain 
3 
= 5, +0 (1), (3.36) 
so that (3.37) 
The equation (3.29) for ‘’, therefore reduces to 
3 2 1 u“ 220 
DEY, = (1+ — (3.38) 


and the particular integral that is not of greater order than unity in the 
Stokes region, and whose terms of this order match Stokes’s solution, is 


= — (3.39) 


As was to be expected, this is the rotational part of Oseen’s solution (2.5). 


3.4. The second term of the Stokes expansion 
When the first approximation to the left-hand side of the Stokes form 
(3.1) of the governing equation is calculated from the leading (Stokes) term 
(3.15), the equation for x, may be written in the form 


Hence we may write, for the time being, 


JAR) = R, (3.41) 
provided we allow for the possibility that the arbitrary constants in the 
integration of (3.40) may be functions of R. Such functions must, of course, 
be of smaller order than 1/R in order that the solution should not contribute 
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to the leading term of the Stokes expansion; and there is no point in 
considering functions of smaller order than unity because these will be 
considered in subsequent terms of the Stokes expansion. 

A particular integral of (3.40) which satisfies the boundary conditions 
atr = 1 andp= +1 is 


3 


‘The general solution satisfying these boundary conditions is then obtained 
by adding to (3.42) the complementary function (3.20). ‘The contribution 
that this general solution makes to ‘¥ is therefore 


R 


+ B,{2R-" +29" +1 — (2n + (2n— 1)R"*3p-"}]0, (3.43) 


Now, the leading term of the Oseen expansion (3.13) has already been 
matched to the Stokes expansion, and it has been shown that F(R) = R 
Hence no term in (3.43) must be of greater order than R. ‘Thus we get 


A, =0, 1 
B,=0 (n>2), } (3.44) 
B, = O()). J 
When p = O(1), (3.43) therefore reduces to 
Ri + 2B, + 0(R), (3.45) 


and B, must be chosen so that (3.45) agrees with the expansion of the 
Oseen term ‘’, at small values of p. From (3.39), this latter expansion is 


so that + 2B, Oy(u) = (3.47) 
‘Thus, since = —du(1—p?), (3.48) 
we obtain B, = —. (3.49) 


The solution for the second term of the Stokes expansion is therefore 


It does not seem to have been emphasized sufficiently in the literature 
that a knowledge of Oseen’s solution (2.5) enables a second approximation 
in the Stokes region to be obtained without a complicated application of 
the whole Oseen technique for obtaining a second uniform approximation. 
Yet the derivation of this second approximation to y% is, in principle, a 


€) 
i 
\ 
=, 
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necessary preliminary for the calculation of a second approximation to any 
property of the flow. The case of the drag coefficient, however, is rather 
special. For the second term of (3.50), being odd in », makes no contribution 
to the pressure-drag or skin-friction on the sphere. But this term is the 
particular integral arising from the inertia forces. Hence Oseen’s theory, 
which neglects inertia forces in the Stokes region, must give a proper second 
approximation to the drag coefficient. The actual calculation of this 
second approximation is particularly simple because the first term of (3.50) 
is a multiple (3) of Stokes’s solution. In this way, we obtain Oseen’s 
result (2.6). 

Another interesting property which may be discussed on the basis of 
the second approximation (3.50) is the formation of ‘eddies’ behind the 
sphere. ‘The first two terms of the Stokes expansien may be written in 
the form 


which, for sufficiently small values of R, vanishes only atr = landu = +1. 
For larger values of R, however, ¢ also vanishes along the (real) curve whose 
polar equation is (according to (3.51)) 


3.52 
and this curve is the boundary of the ‘eddies’. Equation (3.52) shows 
that the minimum value of « occurs at 7 = 1 and is 
nin = R 4° 
so that the ‘ eddies’ first appear at the rear stagnation point and do so when 
R=8. Of course, this Reynolds number is far too large to make estimates 
based on only two terms of the Stokes expansion at all reliable. In fact, 
it cannot seriously be claimed that slow-motion theory gives even a 
qualitative explanation of the phenomenon. A similar conclusion was 
reached by Pearcey & McHugh (1955) in the case of Oseen’s (rather than 
the Navier-Stokes) equation, after a careful numerical evaluation of 
Goldstein’s (1929) exact solution*. 


3.5. Higher approximations 

Since higher-order terms in the Stokes and Oseen expansions are not 
proportional to simple powers ot R, it seems desirable to give a brief account 
of the nature of these terms. This can be done most easily by first 
considering the third term of the Stokes expansion, #5. 


*In this respect, the numerical work of Tomotika & Aoi (1950) appears to be 
seriously in error. 
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we set fx(R) = (3.54) 
and allow the arbitrary constants in the integration for %, to be functions 
ot the Reynolds number, in the manner described for the case of u,, an 
elementary application of the matching conditions shows that the general 
solution for 7, must be of the form 


2 


P 3 
he = —rP+3rlogr— — 


40 
3 logr 7 1 


where C,, are constants of integration. ‘The evaluation of these constants 
proceeds as follows. Both C, and Cy follow immediately from the no-slip 
condition on the sphere. ‘The term in C, makes a contribution 3C, Rp* 
to ‘, so that C, may be evaluated from the small-p expansion of the (known) 
function ‘’,; C; and C, then follow from the no-slip condition. The 
case of Cy, C, and C3, however, is somewhat different owing to the presence 
of the particular integral in r*logr. In fact, the contribution to ‘I’ of the 
term in is 
— Rp? — 3R? log Rp? + 3 R®p log p + C, R’p?)O,(u) + 0(R?), (3.56) 
which contains a term in R?logR. If, therefore, there is no term of this 
form in the Oseen expansion (and it will be shown below that this is indeed 
the case), we must have 
C, =3log R+ O(1). (3.57) 
The constants C, and C, must then be multiples of log R also, in order to 
satisfy the no-slip condition on the sphere. Indeed, it is clear that we 
should replace (3.54) by 
fx(R) = log R (3.58) 
and that the corresponding #, is a finite multiple of Stokes’s solution (3.15). 
Thus, using (3.55) and (3.57), 
‘The proof that there is no term in R?log R in the Oseen expansion is 
very simple. Such a term would have to satisfy Oseen’s equation (3.26) 
and the general relevant solution for D3‘, would be (3.29), in which 
the coefficients 4, are of order unity. But the contributions to D? % would 
then be 


9) 


a 


(Relog Rye!” 4,(4Rr)K, (1) (3.60) 


and (3.32) shows that the order of magnitude of (3.60) in the Stokes region 
is RlogR or greater. Such terms are known not to exist, so all the 4, 


P 
| 
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must be zero. A similar proof applies to the solution of the resulting 
equation DEY 


Hence there is no term in R? log R (or, indeed, in any function of R whose 
order of magnitude lies between R and R?) in the Oseen expansion, and 
(3.58) and (3.59) are a valid representation of the third term of the Stokes 
expansion. 

Nevertheless, the Oseen expansion must ultimately involve terms in 
log R. For the strength of the source-flow outside the wake in the Oseen 
region is known to be proportional to the drag coefficient (Goldstein 1929). 
And, according to the first three terms of the Stokes expansion, the drag 
coefficient is 
R+O(R? 3.61 

(the third term is a new result). ‘The argument shows, in fact, that the 
first occurrence of logR in the Oseen expansion is in the term R®log R. 
‘The non-linearity of the Navier-Stokes equation then shows that both 
expansions must involve powers of log R, and it seems reasonable to suppose 
that both expansions are in powers of R, each term of which is multiplied 
by a polynomial in log R. 

Finally, it is worth noting that, although the technique was not designed 
with such in mind, the Oseen expansion is actually uniformly valid over 
the whole field of flow. Each successive term of the expansion does not 
give a higher uniform approximation, but it is always possible to find 2 
finite number of terms which give, to any required accuracy, a uniform 
approximation to the solution. For instance, a first approximation requires 
all the terms up to, and including those in R’, since, at that stage, the last 
term (the doublet) of Stokes’s solution has been matched. The fact that 
it is possible to choose a finite number of terms in this way provides a 
very simple @ posterior? justification for the formal matching procedure. 


4. GENERAL EXPANSIONS FOR FLOW PAST A CIRCULAR CYLINDER 

In this section, we attempt to obtain higher approximations to the 
velocity distribution for flow past a circular cylinder, and to discover how 
far the method adopted for the sphere in $3 will apply in this case. This 
problem has been considered in outline by Lagerstrom & Cole (1955) 
who introduce Stokes and Oseen variables and obtain Stokes and Oseen 
expansions that follow naturally from the limit processes they adopt. This 
present account considers the problem in rather more detail, in order to 
display more fully the structure of the expansions*. In§4.5 some comments 


* At a late stage in the preparation of this paper, a short paper by Kaplun was read 
at the IX International Congress of Applied Mechanics entitled ‘* Low Reynolds 
Number Flow Past a Circular Cylinder’’. Though by no means a detailed account, 
it appears to present the same conclusions as are reached in the present paper, and 
gives an improved approximation for the drag. 
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are made on the effect of those inertia terms that are transcendentally small 
compared with the expansion parameter. 


4.1. The Stokes and Oseen expansions 


In the Stokes region of the flow, the equation for the stream function 
pecomes 


R 


Vie = -- 4.1 
r e(r, ( ) 
vhere r, 0, and ¢% are as defined in $2.2, and where 
1 
oe ror. (2) 


We assume an expansion of the form 
b = fol 8) + 8) + (4.3) 
where Fr 90 as R+O (4.4) 
and remains bounded at R= 0. As in the three-dimensional case, 
the expansion (4.3) is to satisfy the equation (4.1) and the no-slip condition 
on the cylinder, while the condition at infinity must be replaced by the 


condition that the expansion should match an expansion which is valid in 
the outer region. Again, for similar reasons, we introduce Oseen variables 


p = Rr, VY = Ke, (4.5) 
in terms of which the equation (4.1) becomes 
1 a('¥, V3 ‘F) 


(4.6) 
We now assume an expansion 

+ Fy(R)¥ (p, 6) + Falp, 8) + (4.7) 
where F,.,(R)/F,(R) +0 as (4.8) 


‘The expansion (4.7) is to satisfy the equation (4.6), the uniform stream 
condition at infinity, and, as before, a matching requirement for small 
values of p: that it should match the Stokes expansion (4.3). 


4.2. The leading terms of the expansions 
From the approach described in §2.2, it follows that we may expect 
the leading term J, of the expansion (4.3) to be the Stokes solution (2.7) 
with f)(R) = C, and the leading term ‘Vy of the expansion (4.7) to be the 
uniform stream 
¥, = psin®. (4.9) 


That (4.9) is in fact correct can be deduced in exactly the same way that 
‘fy was obtained for the sphere in $3.2. Similarly, we can say that %» must 
satisfy 
Vig, = 0 (4.10) 
F.M. s 


= 
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since f,(0) is bounded. The most general solution of (4.10) that is anti- 
symmetric about 6 = 0 is 
by = [B,(2r + — 27 + 1-1) ] sin + 


+ > [B,(r?-(n4 + C,(r" — + (n— 1)r—™)]sin 28, 
n=2 
(4.11) 


where B,, C,, (n >1) are constants. When (4.11) is expressed in terms 
of the Oseen variable p the contribution of f)(R)y to Y becomes 


[B,(2p log p — 2p log R— p+ + G(p?R + — 2p) sin 6 + 
> [B,,(p” +2R-n-1 (n + np *R*) 
+C,,(p"R-“ +1 — R"-1 + (n— nth, (4.12 


If this is not to contain any terms of greater order than unity, then 


C, =90, ] 
(4.13) 
C, =0 (n > 2), | 
fy(R) = I/log R. J 
Hence (4.12) becomes) __ 2B, psind + O(1/log R), (4.14) 


and because of the matching condition and (4.9) we get 

B,= —}. (4.15) 
In view of the non-linearity of the Navier-Stokes equations, the term of 
order (log R)! in (4.14) suggests that the functions F(R) will be inverse 


powers of log R. 


4.3. The second term of the Oseen expansion 
Having established that ‘’) represents a uniform stream, we can use 
Oseen’s equation to solve for ¥'}. In terms of the stream function this is 


= 0, (4.16) 
where = pcos@. (4.17) 
A first integral to this equation that is bounded at infinity is given by 
vn, K,(4p)sinnd (4.18) 


n=l 
where X,, are constants and K,,(3p) is a modified Bessel function. Without 
solving for ‘'; we may apply the matching requirement directly to the 
vorticity, V3‘F', in the form (4.18). We know, from the Stokes solution “,, 
that 


tol iby = rlogR sin (4. 19) 


al 


x 
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Writing (4.18) in terms of the Stokes variable 7, we get 


K,(ERr)sin nd. (4.20) 


n=1 
Since the K,,(z) behave as s~” for small values of z, and since RV} must 
match V?:s, we must put 


A, = 1, A, =0 for 1, (4.21) 
and F(R) = (log Ry". (4.22) 
Equation (4.18) now becomes 
= ef? (3p)sin (4.23) 
and this may be integrated (see T’omotika & Aoi 1950) to give 
2 b, (4p) +harmonic function, (4.24) 
where =2K, 1, + (4.25) 


the A,, and J,, being modified Bessel functions. Now ‘tf’, must tend to 
zero for large values of p and so the harmonic function in (4.24) can only 
be of the form 

Y,,p "sin nd, 
‘The matching condition between Ry and ‘’ then requires that the constants 
\,, must vanish for all. Thus (4.24) reduces to the relevant part of Lamb’s 
result (2.13). 


4.4. Higher terms in the Stokes and Oseen expansions 
From what has been shown above, it is clear that 
F(R) = (log (4.26) 
It is not difficult to see that the Stokes solution given by (4.11) and (4.13) 
will not include terms of order (log R) *, (log R)~*, etc. Hence we conclude 
that 
f,(R) = (log Ry" (4.27) 
with the consequent result that 
= for all n. 


‘This means that as far as the Stokes expansion is concerned, the inertial 
terms are never important enough to be considered in the governing 
equations ; their effect is felt only through the outer boundary or matching 
condition. All the arguments applied to the general solution for % will now 
apply to the general solutions for ,, and we can infer that the ¥,, differ 
from one another only by a numerical factor. ‘This means that the Stokes 
expansion has now been reduced to the form, 


b = f(R)bo + OCR), (4.28) 
where 9) = (2rlogr—r+r-')sin# (4.29) 
and f(R) = 4,(log + x,(log R)? =..., (4.30) 
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the x, being constants. So far we have determined z, = —}. ‘The next 
step is therefore to solve tor V,. This is given by the equation 
p @) 


and must vanish at infinity. The calculation is straightforward, the 
solution involving a particular integral from the right-hand side of (4.31) 
and a complementary function given by (4.18). If V3‘, is expressed in 
the form (4.23) and ‘V, in the form (4.24) then we obtain readily 


where g,(R) = O(R 4), } 
g(R) = O(R *), | (4.33) 
= O(R?*) (p> 2). | 
If we write Vo", = eI, then (4.31) and (4.32) give 
= g,,(p)sin (4.34) 
Solving (4.34) by the method of variation of parameters we get 
= k,,(p)sinn@ + complementary function, (4.35) 
u=1 


where k,,(R) is OUR" *) except for n = 1 when it is O(R). By the same 
argument as that used in $4.3, the complementary function in (4.35) must 
be Z, K,(4p)sin @, and hence (4.35) becomes 


V24, = &,(p)sinnd + Z, K,(4p) sin 8. (4.36) 
| 


Finally, on integrating again, we get 
F, Z (3p) nO (4.37) 


where = (4.38) 


and in particular /,(R) = L,(R)+O(R?) where L, is a constant. When 
F(R), is written in terms of the Stokes coordinate +, it is found to give 
a contribution — Z,(log R)'r to f,(R)by. Thus, from the matching require- 
ment, we deduce that 


= ($}—y+log 4). (4.39) 


Next, the vorticity term V3, given by (4.36) may be matched to the 
corresponding vorticity term RV? >, in the way that (4.19) was matched 
to (4.20), to give 

4 = —4(}-y+log4), (4.40) 


4%, being the coefficient of (log R)-? in (4.30). 
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It is not difficult to see how the process outlined above to give ‘’, and 
thence x, may be used successively to obtain all the 4’), z,,.. Each of the Y’, 
will consist of a complementary function 


psin 


R, > ¢,,(3p) 


where R,, is some constant, and a particular integral 
> m,,,(p)sin 


n=1 


where the m,,,,(p) display the same behaviour as the /,,(p) in (4.37). 


4.5. Further effects of the intertia terms 

The expansions (4.3) and (4.7) which have been considered above 
ditfer very greatly from the corresponding expansions (3.5) and (3.13) 
obtained for the sphere. In fact, when (4.3) is written in the form (4.28) 
it is seen that the entire Stokes expansion obtained for the cylinder corre- 
sponds to just the first term of the Stokes expansion for the sphere. Close 
to the cylinder, virtually no account has been taken of the inertia terms, 
while in the Oseen region terms of order R have been neglected. 

From a physical point of view the expansions that we have obtained 
cannot be expected to provide much information. Many important 
characteristics of slow motion will be caused by just those terms of order R 
that have been neglected because R is transcendentally small compared 
with any power of (log R)'. We can attempt, in a formal way, to take 
account of these terms of order R (and also of higher powers of R) by writing 
our expansions for the stream function in the form 


= (log 8) + (log RY 8) + (log RY 8) + 
+ R(log 8) + (log RY 8) + (log RY 8) + 


+ R"(log 8) + (log R) 8) + (log o(1, 8) + 
(4.41) 


for the Stokes expansion, where fy, §;,...8, are integers, with a similar 
form for the Oseen expansions. 

‘These expansions would have to satisfy the same equations, the same 
boundary conditions and the same matching condition as in $4.1. Formally, 
the derivation of the terms y,, and ‘l’y,, would be carried out exactly as in 
$4.2, $4.3 and § 4.4. It is seen readily that the %,, substituted into the 
right-hand side of (4.1) would yield inhomogeneous equations for the y,,,. 
The matching conditions with the boundary conditions would then enable 
the arbitrary constants in the integration to be determined. By a process 
very much like that used for the #),, and ‘,, all the ¥,, and ‘¥’,, could be 
derived successively. It will be observed that only a finite number of the 
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;, and ‘,, would be needed in order to obtain any particular one of them; 
in this sense, the formal process suggested above for determining the various 
functions in (4.41) is a possible one. 

Asymptotic expansions for low Reynolds number flow are used 
extensively to give numerical approximations to certain constants of the 
physical flow. For R as small as 10-*, R is larger than (log R)-*, and so 
from a numerical point of view the term %4, in (4.41) might well be more 
important than %3. ‘Though it is difficult to justify in any analytic sense 
the use of (4.41), it is nevertheless clear that many properties of slow 
streaming motion past a circular cylinder which do not form any part 
of the asymptotic expansions (4.3) and (4.7) might be predicted by using 
a suitable number of terms from the conjectural expansion (4.41). 
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The inhomogeneity of grid turbulence 


By H. L. GRANT and I. C. T. NISBET 
Cavendish Laboratory, University of Cambridge 


(Received 24 December 1956) 


SUMMARY 

‘The intensity of turbulent motion in the wake of square 
mesh grids of bars has been examined experimentally. The 
measurements show that there are considerable departures from 
homogeneity over surfaces normal to the mean stream. Large 
variations of intensity have been observed at a distance of 80 
meshlengths from the grid, where they are diminishing slowly, 
if at all. This inhomogeneity is believed to be the cause of the 
discrepancies between various sets of measurements of energy 
decay in the ‘initial period’. ‘There is also a possibility of error 
in other measurements made in grid turbulence and analysed 
on the assumption of homogeneity. Observations with a grid 
of 2 in. mesh show that z? is significantly less than 12, which is 
consistent with results from other laboratories. 


INTRODUCTION 

‘The simplest type of three-dimensional turbulence that can be imagined 
is isotropic homogeneous turbulence, in which the statistical properties of 
the motion are unchanged by rotation or reflection of the coordinate axes. 
These restrictions have permitted much greater theoretical progress than 
has been possible with more complicated forms (Batchelor 1953). The 
closest practical approximation to this appears to be the turbulence in the 
wake of a square-mesh grid of bars; this suffers from the disadvantage that 
the intensity of the turbulence decays with time and hence with distance 
from the grid, but this departure from spatial homogeneity in the direction 
of the mean flow is found to be small over lengths comparable to the linear 
scale of the turbulence. Since ‘Taylor (1935) first suggested the use of 
grids for generating turbulence, it has been common practice to assume 
that this ‘grid turbulence’ is homogeneous and isotropic in this approxi- 
mate sense, and it has been used almost exclusively in the experimental 
studies which support the theory of isotropic turbulence. 

It had recently become clear from unpublished work by Corrsin at 
Johns Hopkins University and by Wyatt at Manchester that such turbulence 
is not accurately isotropic, and we therefore attempted to determine the 
ratio of the intensities of the down-stream and cross-stream components of 
turbulent velocity in the small wind tunnel in the Cavendish Laboratory. 
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By measuring this ratio behind several grids and over a wide range of 
turbulent intensities, it was hoped to establish standards for the calibration 
of hot wire anemometers. 

During the investigation it was found that the turbulence was not only 
anisotropic but also inhomogeneous. ‘This lack of homogeneity has 
previously been recognized by workers using very fine grids of mesh length 
4 in. or less (Batchelor & ‘Townsend 1948; Batchelor & Stewart 1950), 
but it does not seem to have attracted attention in connection with the 
larger grids. We therefore turned our attention to an examination of 
the degree of inhomogeneity of the turbulence, this being an even more 
serious departure from the assumed conditions than lack of isotropy. 


EXPERIMENTAL ARRANGEMENT 

‘The measurements were made in a wind tunnel with a working section 
15 in. square and 8 ft. long. ‘lhe turbulence was produced by passing the 
air through biplane grids of circular cylinders of mesh to diameter ratio 
5-33. Each grid is made from two plane arrays of parallel, equally spaced, 
circular cylinders. ‘The two arrays touch each other and have the directions 
of their cylinder axes at right angles to each other. Five grids have been 
used and it is convenient to designate them by their mesh lengths, which 
are }in., }in., lin., and (in two cases) 2in. ‘They are geometrically 
similar, except that the intersections of the } in. grid are soldered, and the 
mesh to diameter ratio of the 1 in. grid is slightly greater than 5-33. ‘The 
wires of the } in. grid are under tension. One of the two 2 in. grids, which 
was kindly lent to us by Professor 5S. Corrsin of the Johns Hopkins 
University, is made of wood with small nails at the intersections; the other 
2 in. grid and the 1 in. grid are made of brass rods with no mechanical 
connection at the intersections. 

The mean stream velocity will be designated by U, the intensity of the 
component of turbulent velocity in this direction by v2, and the intensities 
of the turbulent velocity components parallel to the two sets of bars in the 
grid by 72 and w. Wind speeds were either 6-3 or 11 m/sec. 

The hot wire anemometers were of platinum-cored wollaston wire with 
a diameter of 0-00025 cm after etching. A single wire normal to the mean 
stream was used for most determinations of v2 but, in cases where two 
components were measured, an x-wire was used for both. The wire length 
was usually about 0-05 cm and never exceeded 0-1 cm for either single 
or x-wires. ‘The wires were placed ina Wheatstone bridge and the unbalance 
amplified with appropriate compensation for the lag caused by the thermal 
capacity of the wire. The output was fed to a vacuum thermojunction 
supplying a millivoltmeter. In the case of a single hot wire normal to the 
mean stream, this arrangement yields Cu’, where the calibration ‘constant’ C 
is a function of the dimensions of the wire, the mean velocity, the heating 
current, the gain of the amplifier, and the amount and distribution of dirt 
on the wire. To obtain the intensity relative to that at some reference point 
in the flow (or in some other flow) simply requires a reading of the mean 
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square output at the reference point enough often to determine any drift 
in wire sensitivity. ‘This is what has been done in cases where the data are 
given in ‘arbitrary units’; the probable error is then about 1°,. 


INHOMOGENEITY OF THE DOWNSTREAM COMPONENT 

Figure 1 shows the spatial distribution of intensity of the downstream 
component of velocity near the axis of the tunnel behind the $ in. grid. 
The diagram represents a plane normal to the mean flow, and we have drawn 
lines of constant u? covering an area approximately 6 in. square, 80 mesh- 
lengths downstream from the grid. The extreme variation of intensity 
in this example is about 30°,, of the mean, and the maximum gradient is 
about 20°, in }in. No cross-stream variations of mean velocity were 
found with a pitot tube and manometer, which was capable of detecting 
differences of 
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Figure 1. Distribution of 42 (in arbitrary units) over a plane normal to the mean 
flow and at 80 mesh-lengths downstream from the } in. grid. ‘The contours 
are derived from observations at intervals of } in. 


In the region of the grid corresponding to figure 1, the diameter of the 
bars is very uniform and the spacing does not vary by more than 2°,. The 
intersections are soldered but the solder is in no case visible in the projection 
seen by the wind, although the shape of the fillets varies slightly. There is 
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no significant correlation between the intensities in figure 1 and the area 
of the grid openings; nevertheless, the pattern is a property of the grid, 
since when the grid is displaced sideways the pattern follows it, at least 
for small displacements. 
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Figure 2. Variation of (in arbitrary units) along a line normal to the mean flow 
with various grids. 

Similar effects are found behind all the grids which we have examined. 
Figure 2 shows several examples of the variation of the intensity along lines 
normal to the mean stream. ‘The larger meshes produce more uniform 
turbulence, although variations of as much as 6°,, are found behind the 
2 in. grids. 

These patterns are found to change only slowly with distance from the 
grid, although none has been traced closer than about thirty mesh-lengths 
trom the grid. ‘That such patterns, once formed, would be expected to 
persist for a long distance from the grid may be shown by the following 
rough calculation. If an excessively high intensity is being produced at 
a point in the plane of the grid, the subsequent spread of this high intensity 
in the lateral direction can be regarded as a diffusion problem, with the 
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ditfusion coefficient K given by the product of a representative length 
scale L, and a representative velocity (u?)! of the turbulence. If these are 
constant, the spread of turbulent energy in the lateral plane is equivalent 
to heat conduction from a line source. The distribution of excess intensity 


will be like 


— +Kt 
where ¢ is time measured from the instant the fluid left the grid and y is 
distance normal to the mean streamline which passes through the source 
of high intensity. The half-width of the distribution is 2(K?)?, which is 
equal to 2{Lx(u?)!/U}. If we assume that L/M = 3, and (u?)!/U = 0-014 
at x/M = 50, where M is the meshlength, then 2(At)! = 0-7. If such 
sources are at least a meshlength apart, the turbulence has not had time, 
at 50 meshlengths from the grid, to become homogeneous through adjacent 
‘wakes’ overlapping. It may also be noted that until adjacent ‘wakes’ 
do overlap appreciably, the quantity 1/At in the equation is decreasing as 
v3, so that the excess intensity is decreasing rather slowly along a 
streamline. 


INHOMOGENEITY OF THE CROSS-STREAM COMPONENTS 


The intensities of the cross-stream components of velocity have similar 
variations, but their spatial distributions are not obviously correlated with 
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Figure 3. he distribution of the intensities of the three components of turbulent 
velocity near the point X in figure 1. The letters on the abscissa refer to 
points located with respect to X as shown in the inset. The intensity of each 
component is taken as 100 at _X. 
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each other or with that of the intensity of the downstream component. 
An example from the 3} in. grid is given in figure 3, which shows the 
distribution of intensity about the point X in figure 1. This point was in 
the region where wv? was most uniform, yet quite large variations were found 
in the other components. The obvious anticorrelation between v? and a? 
is probably not typical, but we have not measured all three components 
in enough cases to be sure. We do have a considerable number of cases 
where the intensities of the downstream and one cross-stream component 
have been recorded, and there is usually no obvious correlation between 
the two curves. 


ABSOLUTE DETERMINATION OF INTENSITY 

It is clear from these observations that the ratios of the intensities of 
the ditferent velocity components are not constant, so that the turbulence 
cannot be isotropic. We have made a series of determinations of u2/U?, 
z?/U2 and v/v”, at 30 meshlengths downstream from one of the 2 in. grids 
at a point which will be referred to as reference point A. 

For the determination of u?/U?, a single hot wire normal to the stream 
was calibrated by varying the wind speed in the absence of a grid and 
observing the variation of heating current required to keep the wire at 
constant temperature. From this, v?/U? was obtained indirectly by using 
a single hot wire which could be operated in three positions, perpendicular 
to the mean stream and at 45° on each side of it. ‘This arrangement gives 
signals proportional to the mean squares of u/U, (u+v)/U and (u—v)/U, 
the ratios of the sensitivities being determined by the mean heating currents 


U (cm’sec) | Component Value Remarks 
630 u?/U? 45 «104 Measured. 
630 0-78 +0-02 Measured. 
630 v2/U? 3°5 xdo- Deduced from above figures. 
1100 ne/U? 43 x10-! Measured. 
1100 0-72 +0-01 Measured. 
1100 Deduced from above figures. 
1100 0-105 0-005 Measured values of 
varied between —0-08 and 
—0-12 over this lateral plane. 


Table 1. ‘Turbulent intensities at reference point 4. 


in the three positions. Accurate alignment with the mean stream can be 
etfected by equalizing the mean heating currents in the two 45° positions, 
and hence values can be obtained for w2/z2 and uz/u®?. An attempt was also 
made to determine v/U? directly, using an x-wire which was calibrated by 
rotating it in its own plane about an axis normal to a steady flow, thereby 
introducing a known v-component to the stream past the wire. This 
method did not give consistent results and was discarded. 
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In view of the inhomogeneity, the results (see table 1) are useful for 
calibration purposes only in this particular tunnel, but the value of v?/u? 
is of more general interest because it is smaller than unity by about 25°. 
This is consistent with the values found by Corrsin and Wyatt, and is 
significantly different from unity since the maximum variation of uw? and v? 
is about 3°, from the mean in the case of the 2 in. grids. 


DECAY OF TURBULENT ENERGY 
The decay of turbulent energy in the ‘initial period’ has been the 
subject of a number of experimental studies and much theoretical discussion 
{see Batchelor (1953)). We do not know how the rate of decay of energy 
along a streamline of the mean flow depends upon the lateral intensity 
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Figure 4. Decay curves taken at different lateral positions in the tunnel with the 
1in. grid. No great care was taken to see that the hot wire followed a mean 
streamline. 


gradient, but it is certain that inconsistent results can easily be obtained 
through failure to make the intensity determinations for a decay curve 
along a single mean streamline. As an illustration of what can happen, 
we present in figure 4 two determinations of U2/u? as a function of x/M 
for a wind speed of 6-3 m/sec. The same grid was used, but the hot wire 
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was moved down the tunnel at two different lateral positions without taking 
particular care to see that it followed a mean streamline. As long as the 
hot wire follows the same track, each of these curves is quite reproducible. 


OTHER GRIDS 

Brief surveys were made of the turbulence behind two other types of 
grid. One consisted of a sheet of metal 0-02 in. thick in which 3 in. holes 
had been punched at regular intervals in a hexagonal pattern leaving a 
blockage ratio of 0-48. The edges of the holes were very neat and the 
diameter did not vary by more than 0-001 in. ‘The distance between the 
holes was uniform to within 0-4”, in one direction but there were variations 
of +°,, from the mean in the perpendicular direction. ‘he intensity 
behind this grid was very irregular in both directions but slightly more so 
in the direction along which the hole spacing was not uniform. At one 
place the intensity changed by a factor of two in 2in. ‘The other grid was 
a single array of parallel rectangular prisms of section 3-7 cm x 0-4 cm, 
spaced 19cm apart (Stewart & ‘Townsend 1951). Here the variations 
over mosi of the tunnel section were within only 2 or 3°, although at one 
point the intensity rose to 15°, above the mean. 


DiscussloN 

‘The possibility of appreciable inhomogeneity has not been considered 
when previous measurements of the properties of grid turbulence have 
been interpreted with the aid of the theory of isotropic turbulence, and it 
is necessary to consider whether or not appreciable errors may have occurred. 
In most measurements one of the independent variables has been distance 
downstream, and it can now be seen that unless care is taken to prevent 
it, the probe may wander from a mean streamline on which the turbulence 
perhaps has a high value of uv? to a mean streamline on which u? is lower 
than the average value over a lateral plane. We believe that this is the chief 
cause of the unreliability of our decay curves, but it is also possible that 
there are real ditferences between the decay curves at different parts of the 
flow and that these contribute to the variations among the observations. 
It seems likely that these differences would be small where the lateral 
intensity gradient is not large. If this is so, then a significant determination 
of the rate of decay could be obtained by following a mean streamline. 
‘To follow a mean streamline with sufficient accuracy would not be difficult 
with the 2 in. and larger grids; with which the scale of the variations is 
correspondingly large; however, a very long wind tunnel is needed to give 
access to the interesting final period of decay with these grids. 

Studies of the ‘approach to isotropy’ appear to be very unreliable if 
made in this kind of turbulence, as we have obtained curves of v?/u? against 
distance from the grid which sometimes had a positive slope and sometimes 
a negative slope. Most such studies, however, have been conducted with 
turbulence which has deliberately been made highly anisotropic (Townsend 
1950, 1954; Batchelor & Stewart 1950), and it is hard to estimate the 
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effect of an initial inhomogeneity in the ratio v2/u?. In these cases any 
error is more likely to have arisen from the assumption, used in calibration of 
the hot wires, that o? = w? in grid turbulence in the initial period of decay. 

Similar comments could be made_about functions of the velocity 
derivatives. In the one case in which (du/dx)* was measured, it was found 
to vary in a similar manner to uz, but it is not safe to conclude from this 
one observation that the spatial variations of the two quantities are always 
correlated. 

Of the correlation functions, f(r) = R,,(r,0,0)/u? is the one of most 
interest in isotropic turbulence. It is probably not much affected by spatial 
variations in the intensity unless the proportional distribution of energy 
among the various scales of turbulence also varies in the lateral direction. 
We have no information about this. Another correlation component which 
has often been measured is g(r), one form of which is R,,(0,7,0)/u?. The 
most closely related quantity which can be defined in a flow which is not 
homogeneous is the correlation coefficient R,,(0,7,0)/(u?w’?)!, which in 
fact is what has actually been measured in most cases. (u2 and uw are the 
intensities at the positions of the two probes.) This is determined without 
obtaining the two intensities explicitly, and is independent of variations of 
intensity although it would be affected by variations of scale. 

There are certain measurable properties of the turbulence which are 
likely to depend only on the gross geometry of the grid and which could 
be expected to hold for wind tunnels and grids other than those with which 
the measurements were made. One such quantity might be the mean, 
over many points in a lateral plane at a fixed distance from the grid, of the 
local time average of wu. ‘This combined space and time average might be 
an acceptable approximation to the intensity in the theoretical homogeneous 
turbulence. However, this mean is derived from a velocity distribution 
which may not be independent of the detailed initial conditions prevailing 
at the region of formation of the turbulence (detailed geometry of the 
intersections, level of turbulence in the incident stream, etc.), so higher 
moments might be very different from those of homogeneous turbulence 
even if defined in the above way. 

For a proper understanding of the relation between real grid turbulence 
and truly homogeneous turbulence, it would be necessary at least to make 
a study of the decay along mean streamlines. Even more useful would be 
to find a way of making turbulence which is really homogeneous if not 
isotropic. 

CAUSE OF THE INHOMOGENEITY 


Very little can be said about the cause of the phenomenon. It may be 
that the point of separation on the cylinders forming the grid is sensitive 
to very small changes in the geometry of the intersections, or to surface 
roughness, although this is not likely to be a factor with the grid made by 
punching holes in a sheet. One of the authors is at present conducting a 
study of the flow in the region immediately behind the grid and it is hoped 
that this may throw some light on the subject. 
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SUMMARY 


The partial differential equations which describe steady flow 
of fluid in saturated homogeneous permeable solid material under 
non-isothermal conditions are stated. From these are derived the 
equations for flow of liquid (in particular, water) using suitable 
approximations and making use of empirical laws when necessary. 

It is then postulated that the only ‘ ponderomotive’ (i.e. mass- 
moving) forces present are those due to thermal expansion effects. 
Free convection results. An approximate solution of the equations 
is attempted for plane flow by means of classical perturbation 
methods, the temperature and stream-function variables being 
represented by power series in a convection parameter proportional 
to the Rayleigh number. 

A numerical example of the method, with boundary conditions 
based on a geothermal area at Wairakei, New Zealand, is given. 
The results show features which are in fair agreement with 
temperature measurements made in the area, and it appears that 
the convection parameter 7 is of the order of 10. 


1. INTRODUCTION 


In some situations associated with geothermal activity, it is possible 
that the flow of ground water of meteoric origin is influenced by convection 
currents due to differential heating. Heat is transported both by the 
convecting liquid and by thermal conduction through the saturated permeable 
earth. When an area of thermal activity is under investigation for its 
potentialities as a power source or for other reasons, the relation between 
ground-water movement and temperature distribution is of practical value, 
since the latter can be measured relatively easily by means of exploratory 
boreholes. 

The aim of the present study is to deduce approximate solutions for 
the flow field and temperature distribution when heat conduction and free 
convection of water in the permeable material are important. As would 
be expected, the problem is closely related to the problem of free convection 
in a viscous fluid (Goldstein 1938; Batchelor 1954), in which the apprepriate 
convection parameter is the Rayleigh number, which, in the liquid case, is 

Re a( — (1) 
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Here 7,, 7, are two reference temperatures on the absolute scale, « is the 
coefficient of linear thermal volume expansion of the liquid, g is the 
acceleration of gravity, d is a representative linear dimension, « is the thermal 
ditfusivity of the liquid at temperature 7), and vp is its kinematic viscosity 
at 

When the liquid flows slowly through permeable material, motion is 
resisted according to the law of Darcy (1856), which states that the hydraulic 
gradient is proportional to the fluid velocity and to its viscosity, and is 
inversely proportional to the permeability. In the equation of motion, the 
Darcy resistance term replaces the Navier-Stokes viscosity term. ‘The 
corresponding convection parameter 7 can be defined by 


R, (2) 


where K,,,/K,, is the ratio of the thermal conductivity of the solid—liquid 
mixture to the conductivity of the liquid, and « is, as before, the thermal 
diffusivity of the liquid. The permeability of the solid medium is k = NA?, 
A being a linear dimension proportional to the particle size, and N a 
dimensionless constant dependent upon the geometrical shape of the 
particles. 

Usually the ratio »/R is very small, as the permeable solid offers 
considerable resistance to the flow; and large temperature gradients are 
necessary before appreciable heat transport by convection is observed. 
This suggests the use of series expansions of the dependent variables 
(temperature and stream function) in powers of the parameter 7, analogous 
to the known technique of expansion in powers of the Rayleigh number 
as used in certain problems of free convection in fluids. The method has 
limited usefulness in the latter case, since the series are convergent only 
for small values of R (Batchelor 1954). However, when the heat-transfer 
process is dominated by conduction, as in the case considered here, the 
expansion method should prove valid for a moderately large range of 
temperature differences. 

It will be assumed that the permeable solid medium is homogeneous 
and isotropic in its physical properties, including fluid permeability and 
thermal conductivity, both of which will be assumed to stay approximately 
constant with changes in temperature. 

Because of the wide temperature range employed, it is necessary to 
adopt certain empirical results in order to describe the density and viscosity 
of the liquid (water). As usual, the effect of pressure changes upon the 
density will be neglected. It remains to describe the effect of thermal 
volume expansion, which follows a non-linear law in the case of water. 
For convenience, the approximate empirical expression will be taken in 
the form of a polynomial relating density p and temperature 7 (°C), viz. 


p = 0-9969{1 — 3-17 10-4(7— 25) — 2-56 x (3) 
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Comparison with values of the density p,,, for saturated water tabulated 
by Dorsey (1940) shows that the density difference (p—0-9969) is given 
to an accuracy of 2% in the range 25°C to 250°C. Fora lesser temperature 
range, a linear expansion law might prove to be adequate. 

Changes of liquid viscosity with pressure changes will be neglected 
(Dorsey 1940). However, the liquid viscosity is a rapidly-varying function 
of temperature which requires to be approximated analytically. The 
exponential expression due to Andrade (1934) suffices for most liquids, but 
it is not very tractable and water does not conform to that law. Consequently, 
for water the approximate empirical expression 


v = 0-332(7 + 13)-'cm? sec“! (4) 


for the kinematic viscosity v, based on tables (Dorsey 1940), will be adopted. 
In the range 15°C to 225°C, the expression (4) is accurate to within 5%. 


2. STEADY-STATE DIFFERENTIAL EQUATIONS 
Let the temperature field be continuous throughout a given region of 
the fluid-saturated homogeneous permeable medium, and let steady-state 
conditions exist, so that all partial time derivatives are zero. 
As usual, the steady-state equation of continuity (Hubbert 1940) is 


V.(pq) = 9, (5) 


where V is the Laplacian operator, p is the fluid density, and q is the flow 
vector expressed in units of fluid volume crossing unit area in unit time. 
(It is to be noted that q is related to the fluid particle velocity vector v 
by q = ev, where « is the porosity of the medium.) 

The linear law of motion due to Darcy (1856) can be written as 


~Vp-8+ = —eq. Vq, (6) 


if p, v are the pressure and kinematic viscosity of the fluid, g is the vector 
of gravitational acceleration, and k is the permeability of the solid medium. 
The right-hand side of (6) consists of an inertial term which is negligible 
for the very low Reynolds numbers considered here. Although this term 
is retained in (6) in order that the conventional derivation of the energy 
equation (7) should be applicable, it will be neglected in subsequent sections 
of the paper. 

When non-isothermal conditions apply, both p and v will be one-valued 
continuous functions of the absolute temperature 7. The assumption will 
be made here that equation (6) still holds under these conditions. 

An equation for the flow of energy can be derived also by following 
closely the method for viscous fluids (Goldstein 1938, ch. 14), with the 
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equation (6) replacing the usual Navier-Stokes equation. It is readily 
shown that the steady-state differential equation of energy transport is 


Jepq NT + pV .q— = (7) 


in which K,,, is the thermal conductivity of the fluid-saturated permeable 
medium (assumed constant hereafter), / is the mechanical equivalent of 
heat, and ¢ is the specific heat per unit mass of the fluid at constant volume. 
Equation (7) differs from the energy equation of a viscous fluid only in the 
form of the viscous dissipation term and in the substitution of the flow 
vector q for the fluid particle velocity vector. 


3. APPROXIMATE EQUATIONS FOR LIQUID FLOW 
Now let the fluid be a liquid, assumed incompressible, whose density 
depends on temperature according to the law 


p = —a(T— Ty) — B(T— 
in terms of the parameter 
= T- 
(9) 
Here py, «, 8 are constants, and 7%, 7, are two representative absolute 
temperatures (suggested by the empirical formula (3) for water). Therefore, 
by the equation of continuity (5), one can define a solenoidal vector 


qo = — — — B(T, — Ty)?6?} 
such that V .(pq) = po V -G = 0. (10) 


When converting to dimensionless variables, it is convenient to take 
a length unit d, and to write V (as before) for the new dimensionless 
Laplacian operator. Also, let 


(11) 
where g’ is a unit gravity vector in the dimensionless system, and let the 
dimensionless variables ¢, H, o be defined by 


= Qo (12) 
(K,,/K,,)VH (VP/po g)R/(xv9), (13) 
Vite, (14) 


where « = K,,/py¢ is the thermal diffusivity of the liquid, and v, is the 
value of v, both at 7’ = 7). From equation (2), 7 is given by 
ko(T, — Ty)gd 


(K,,,/K,,)9 


(15) 


For liquid flow through permeable material under moderately low 
pressures and at low velocities, as in near-surface ground-water movements, 
the work done by compression and viscous dissipation is assumed to be 
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small. It follows that the second and third terms on the left-hand side of 
equation (7) may be neglected. 

With the above definitions and approximations, it is foun. that sub- 
stitution of (8) to (15) into (5) to (7) leads to the following dimensionless 
ditferential equations for liquid flow: 


V.o=0, (16) 
VA + + (B/2)(T) — +06 = 0, (17) 
V0 = V0. (18) 
Equations (16), (17) and (18), together with the relation 
o = o(8), (19) 


constitute a system of six differential equations to determine the six 
unknowns 0, o, H and the three components of ©. 


4, FREE CONVECTION WITH PLANE FLOW 
In this section it will be assumed that the predominant fluid motion is 
due to free thermal convection. ‘Therefore, it is convenient to eliminate 
the pressure terms by taking the curl of (17). Some simplification can be 
introduced by assuming plane flow, for which one can write 


—Vx (jp) (20) 
by virtue of (16), j being a dimensionless unit vector normal to the flow 
plane, and %a scalar stream function. ‘Then, from (18) and (17) respectively, 
the appropriate equations of energy and motion can be written as 


= —[V8, Vu, j], (21) 


= ntl + — V4, (22) 
together with (19). (Use has been made here of the vector notation 
a.bxe = [a,b, c}.) 

From equation (4), it is clear that (19) can be written 
o(9) = {1+a(7,— (23) 
where a = 1/(T,—260) (24) 
when the liquid is water. 
On the boundary, it will be assumed, for the present, that 6 obeys the 
inhomogeneous mixed conditions 
d0/du+C,6 = C,, (x = r°), (25) 


where r is a general position vector, C,(r*), C,(r‘) are real functions of 
the boundary position vector r’, and 0/du signifies differentiation normal 
to the boundary. Further, it will be assumed that ys obeys the homogeneous 
conditions 


ys /du+ Czy = 0, (r = r*), (26) 
where, usually, either 1/C,(r8) = 0 (Dirichlet conditions) or C,(r%) = 0 
(Neumann conditions). Solutions of (21), (22), using (23), with the 
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boundary conditions (25) and (26), must be of such form that, as 7 tends 
to zero, % tends uniformly to zero everywhere, and 6 tends uniformly to 
the steady-state conduction solution. 

Although the stability of the flow patterns represented by these solutions 
is an important question, it is outside the scope of this discussion, which is 
restricted to certain specific solutions in the case of slow steady motion. 

The relative importance of convection and conduction in transporting 
heat is indicated by the magnitude of the parameter 7. When 7 is large 
and the convection component predominates, a distinct boundary layer 
could exist. ‘The usual approximations are then valid. However, a 
treatment along these lines is not needed here, because fairly small values 
of » are expected in most geophysical situations. 

When the latter condition that 7» is small applies, it is feasible to seek 
an approximate solution by means of perturbation expansions, as indicated 
previously. Let the dependent variables 6, % be expanded in the power series 


(27) 
(28) 


for which non-zero radii of convergence exist throughout the region in 
which a solution is required. ‘Then, using (23), it is possible to expand o 
in the power series 


o= > (2) 
t=0 
where = {1+ a(T,— (30) 
and o,/0, is equal to the cofactor of the (1, ¢)th element in the infinite triangular 
= 1 


in which use has been made of the Kronecker delta notation. 

When (27), (28) and (29) are substituted into (21) and (22), and the 
successive coefficients of 1, 7, 77... are equated to zero, the following sets 
of equations result: 


V26, = 0, (32) 
and V4, = — (33) 
s=1 
(a, Vis,,) = a( 7 2 2 + Ig VO, J] + 


n—1 


f=() 


forn > 1. 
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For the boundary conditions, it is evident from the substitution of (27), 


(28) into (25), (26) that 


for all n >0, and that 
C3 0, (r = (36) 


forn >1. 

With the above specifications, one can then define a pair of Green’s 
functions F(rlr,,), G(r|r,,) (with position vectors r, r,) according to the 
differential equations 

= —4rA(r—r,), (37) 


= —4rd(r-r,), (38) 
0 


which make use of the notation of the Dirac delta function on the right-hand 
side. ‘These equations have the homogeneous boundary conditions 
oF /du+C, F = 0, = 2°), (39) 
0G/cu+ = 0, = (40) 
so that the formal solution for 6, is a surface integral, and the formal solutions 


for the 6, ¢,, (w >1) can be written down in terms of integrals over the 
volume V of saturated permeable material. ‘These integrals are 


and = — Flright-hand side of (33)} dV, (42) 
JV 
= ~ | Glright-hand side of (34)} dV, (43) 
JV 


forn>1. Here dS,, dV,, represent the surface and volume integration 
elements respectively. Comparison with (33) shows that the right-hand 
side of (42) involves solutions up to @,_,, %,,, and comparison with (34) 
shows that the right-hand side of (43) involves coefficients up to 4,,_;, %,_1- 
Hence it becomes possible in principle to solve the linear equations (32) 
to (34), or the corresponding integral equations (41) to (43), in appropriate 
successive order for 95, 45..., each in terms of the preceding 
coefficients. 


5. CALCULATION OF PERTURBATION COEFFICIENTS 

The practicability of the above perturbation scheme is now tested by 

means of a numerical example. A two-dimensional case is considered, 

and iterative (relaxation-method) solutions are obtained for the differential 

equations (32) to (34) expressed in finite-difference form. Fairly coarse 

nets are employed, and only the first few coefficients of the series, up to 
iby, 4, are calculated. 

The choice of boundary conditions for the model has been influenced by 

a geological situation which appears to exist at one of the geothermally-active 
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areas of Wairakei, New Zealand. In figure 1, EG represents the ground 
surface and HAFG, DAB represent sections of a sheet of ignimbrite 
(an igneous formation) faulted along the line BF. Non-fractured ignimbrite 
has the property of very low permeability, so that fluid flow may be ignored 
in this material. Above the ignimbrite, and bounded by EFAD, are 
formations which possess appreciable fluid permeability. ‘These formations 
appear to be saturated with meteoric water which, although at very high 
temperature in the deeper regions, is held in liquid form by the high pressure. 
Overlying this, and close to the surface EG, is a layer of mudstone which, 
in the calculations here, has been assumed to be impermeable. The region 
below BAH is believed to be a reservoir of trapped stream or superheated 
water at a temperature of about 250°C, and some of this fluid may escape 
through the fault fissure at A into the upper permeable layer. The permeable 
region above AD has been drilled to a depth of 1000 m without encountering 
ignimbrite, so that the existence of the ignimbrite continuation DAB is not 
definitely established. While the angle of dip of the fault BF is uncertain, 
it probably lies in the range 75-90, and has been taken to be 90° here. 


G E 
MUDSTONE 
| 
IGNIMBRITE d 
PERMEABLE | 
FORMATIONS | 
H D 
7-250 C A 
IGNIMBRITE 


Figure 1. Hypothetical geological structure based upon a geological structure 
believed to exist at Wairakei, New Zealand. 


Along BAH and EG, the temperatures are taken to be 250° C and 25°C 
respectively, giving 7) = 25+273 = 298° A, and 7,—7) = 225°A. The 
length unit d is taken as equal to the spacing between AD and EF. For 
the scale of this example, d = O(10°). Other values for parameters are 
obtained from equation (24), which gives a= 1/38( A)“, and from a 
comparison of equations (3) and (8), which shows that x = 3-17 x 10-4(" A) |, 
B = 2:56x 10-% A)-*. The boundary conditions are then as follows. 


On EG: 0=0; 


on HA, AB: 0=1; 

on BC, CE, GH: 0@/cu = 0, where d/éu signifies the derivative normal 
to the boundary ; 

on AD, DE. EF, FA: = 0. 
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The condition c6/eu = 0 on BC is only an approximation to the true boundary 
condition. Also, it will be noted that fluid flow into the region ADEF at A 
has been neglected. 

For the region ADEF, the equations (21) and (22), when expressed in 
two-dimensional form, apply. However, for the regions ABCD, HAFG, 
% can be assumed to be identically zero, and 6 obeys Laplace’s equation. 

In the absence of more precise information, the thermal conductivity 
has been assumed constant over the entire plane, so that at the boundaries 
FA, AD between differing media there is no discontinuity in 9 or in its first 
normal derivative. As there are second-derivative discontinuities in @, the 
isotherms show changes in curvature on the boundaries FA, AD. 

Figure 2 gives the numerical values obtained for the finite-difference 
solutions for 6), 0,...0,, and figure 3 gives the values for y,...74,. In the 
calculation of these solutions, several approximations were involved, as 
follows. 

(a) A boundary singularity in 6 exists at the point A, and partial com- 
pensation for its effect was introduced (Woods 1953) during the calculation 
of 4. As the modification was not large, the effect of the singularity was 
ignored during the calculation of higher 6-coefficients. 


E 
= T T =| 
= 0 | 
10°, = 

-39 49 45 26 -12 =3] 
31 3 21 -7 i 4 4 3 2 
19 -14 5 5 4 ° 
-§ 2 -2 -2 | 
| 
-119 -73' -34 -22 -13! -6 
-87 -79 - 35 -5 9 12 3 
-24 7 15 12 3 (6) 
-8 14 12 1 -4 -5 -4 -2 “1 
- lel - 156 -52 -21 12 -5 
117 - 74 - i ifs} 1e 2 8 4 
-53 8 30 24 13 2 
20 438 2 -7 -@ -4 -2 -1 


Figure 3. Numerical values obtained for the coefficients %,, #2, 3, y4 respectively, 
tabulated for each point of the relaxation net. 


(6) On the boundaries FA and AD, the equations obeyed by the 
coefficients ,, (n > 1) undergo a transition from (33) to Laplace’s equation. 
The treatment of this type of singularity in a finite-difference solution 
involves taking half the contribution due to the right-hand side of (33) at 
mesh points lying on the discontinuity, and in refining the net close to the 
discontinuity. Since the refinement of the net has been omitted for this 
example, errors of approximation must be expected in the vicinity of the 
boundaries FA, AD. 

(c) The calculation of 4, and %, from (34) is approximate in that the 
terms under the first double-summation sign (which are of order m) on the 
right-hand side have been neglected in comparison with the other right-hand 
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terms (which are of order n—1). The resultant errors introduced are 
negligible in the regions of highest temperature, rising to a few percent 
in the cold regions where o is large and is small. 

Owing to the extremely high temperature differences applying on the 
boundaries, the magnitude of the coefficients does not decrease very rapidly, 
and it is necessary to consider the number of terms required to give reasonable 
accuracy, for a given value of the parameter 7. From figures 2 and 3, 
there are indications that the sequences (6,,) and (%,,) are oscillatory, whence 
a reasonably satisfactory estimate should be obtainable by choosing a 
suitable point to truncate each series, e.g. after a change in sign of the last 
coefficient. 

Figure 4 is plotted from the finite-difference solution for 9, and corre- 
sponds to the temperature field due to conduction alone. 


mir 


Figure 4. Field of the temperature parameter 65, corresponding to the conduction 
field in the geothermal model. 


In figure 5, a value of 7 = 7 has been chosen to illustrate the alteration 
in the 9-field with convection present, the series being taken as far as the 
term in 7°. (The effect of adding the term in » is shown by the broken 
lines.) Three departures from the conduction case may be noted: (1) a 
crowding of the isotherms towards point A from the right, corresponding 
to a decrease of temperature in the vicinity of AD; (2) an increase of 
temperature in the vicinity of FA; and (3) a marked horizontal temperature 
gradient as boundary AD is approached from above. If » were increased, 
these phenomena would be intensified, and the isotherms would tend to 
form a ‘mushroom’ situated to the right of the boundary FA. 
Phenomena (1) and (2) have been observed in the data obtained from 
temperature measurements in boreholes at the Wairakei geothermal area, 
and comparison of the temperature gradients with theoretical results 
indicates that 7 = O(10). However, (3) has not been observed, so that, 
although the presence of a convection process appears highly probable, 
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there is no confirmation of the existence of the ignimbrite sheet DAB. 
‘There are indications, also, that the mudstone layer EG has appreciable 
permeability, perhaps due in part to erosion by surface streams. 


G 


Figure 5. Field of the approximate temperature parameter @ defined by 
= 05+0, n?--0, 7°, 
for 7 = 7. The broken-line isotherms indicate the modifications to the 
field when a further term (0, *) is added to the series. 


™m 


A D 
Figure 6. Field of the approximate stream function % defined by 
b= 7? 
for 7 = 7. The broken-line isopleths indicate the modifications to the field 
when a further term (4 7‘) is added to the series. 


In figure 6 are plotted values of the series for the stream function ;/, 
up to the 7° term (full lines) and the y+ term (broken lines), for » = 7. 
Modification by the fourth-order term is in the direction of decreasing +, 
which follows from the sign change in the higher coefficients. Jt will be 
noted that the region of main circulation is concentrated near the cu .er A, 
where a large temperature gradient exists, and where the high temperature 
results in a low fluid viscosity. 
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An estimate of the heat output O at the surface EG (per unit width of 
boundary parallel to the axis of convection) for the boundary conditions 
given is expressible in dimensionless form in terms of the Nusselt number VN: 

O 
|, 
where d/ represents an element of length of the boundary EG. Clearly, 
N is a function of y, and it is found that 
N(n)/N(O) = 1-—0-52 x + 4-25 x 10-4? 5-05 x 10-593 + 2-76 x 
(45) 
N(O) being the Nusselt number for conduction alone, with the given 
boundary temperatures. When » = 7, N(n)/N(O) = 1-04. 

For the boundary conditions of this example, it is evident that the 
chosen value of 7 for the parameter y tests the perturbation method 
practically to the limit of its usefulness. A low maximum value of 7 is to 
be expected, since the temperature difference (7, — Ty = 225° A) is so large. 
Near 250°C, the density of saturated water is about 0-8 of the value at 
25°C, and the kinematic viscosity is less than } of the value at 25°C. 
Nevertheless, these extreme conditions are known to exist in at least one 
geothermal area, and indications are that the value of » appropriate to the 
area does not greatly exceed the value chosen in this example. 
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SUMMARY 


A first-order relationship between changes in area and shock 
strength is derived for the case of a shock moving through a small 
area change in a channel. By integration of this relationship the 
area of the channel is obtained as a function of the shock strength 
in closed form. ‘This result is interpreted as giving the average 
strength of a shock at a given time as it moves along a channel 
of arbitrary shape. 

By suitable choices of the shape of the channel, descriptions 
of converging cylindrical and spherical shocks are obtained. 
These descriptions are found to be in close agreement with the 
similarity solutions valid near the points of collapse of the shocks. 
The reason for such good agreement is examined. 


1. INTRODUCTION 


The motion of a converging cylindrical or spherical shock wave in a 
perfect gas with constant ratio of specific heats has been studied by Guderley 
(1942) and Butler (1954). ‘The shock moves into a uniform medium initially 
at rest. As the shock converges it becomes stronger, until it collapses on 
its axis or point of symmetry, where its strength is singular. Thus in the 
neighbourhood of this axis or point of symmetry, which is taken as the 
origin, the ‘strong’ shock conditions are applicable. In these conditions 
the pressure in front of the shock is neglected in comparison with the 
pressure behind the shock. Guderley’s solution uses these and therefore 
is applicable only in this neighbourhood. ‘This simplification leads to a 
boundary condition for the flow behind the shock wave which permits a 
similarity solution to the problem, including a description of the resulting 
outward-going shock. ‘The similarity variable is R'%/t, R being distance 
from the origin, ¢ the time measured from the instant at which the shock 
reaches the origin, and x a constant depending on y and j; where y is the 
ratio of specific heats of the gas and j has the value 1 for the cylindrical 
shock and 2 for the spherical shock. In this solution the pressure ratio 
across the shock is proportional to R?*°-)*, For y= 5, j = 1, Guderley 
calculates 2(x—1)/x as —0-396, while Butler gives the result correct to 
six figures for y = %, 4, ? and j = 1, 2 with the aid of an electronic computer. 
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An alternative approximate solution to the problem is given in this 
paper. A description of the motion of the shock is given which is not 
restricted to the neighbourhood of the origin, as is the work already 
described. ‘The treatment is based on the description, due to Chester 
_ (1953, 1954), of a shock passing along a channel consisting of two uniform 
sections of nearly equal cross-sectional area joined by a section of varying 
area. Section 2 contains a simple derivation of the asymptotic solution 
after a large time which agrees with the asymptotic form of Chester’s 
solution. After linearization of the problem with respect to 5A, the small 
area change in the channel, the ratio 5z/5A is found in terms of A and g, 
dz being the change at large time of the shock pressure ratio z due to the area 
changeédA. After taking the limit 6A — 0 and integrating, a relation between 
A and z is found in closed form. ‘The relation gives an approximate 
description of the motion of the shock in terms of the area of the channel, 
which may be used for finite continuous area changes. ‘Iwo assumptions 
have been made in deriving this relation. The first is that the effect of 
reflected disturbances generated by the shock may be neglected. In §3 
this assumption is examined for the particular cases of converging cylindrical 
and spherical shocks near the origin, where an exact solution is available 
for comparison. The motion of the shock is found to be only very slightly 
modified by consideration of these reflected disturbances. The work of 
R. B. Payne, described below, suggests that the error due to the neglect of 
these disturbances is also small for symmetrical shocks away from the 
origin. ‘The second assumption is that the steady state solution of Chester’s 
problem, valid only for large time, may be used. ‘This assumption is given 
weight by the work of Chester. He shows that although the shock is not 
uniform across the channel after entering the area change, the change in 
shock strength averaged at any one time over its area is proportional to the 
change in area of the channel. ‘Thus when the area change is passed, the 
average shock strength is constant, there being only local changes in its 
strength tending to make the shock uniform. Hence the integration of the 
differential equation derived from an asymptotic solution will give an average 
shock strength in terms of the area. However, for the special symmetrical 
flows of converging cylindrical and spherical shocks, the shock strength 
at any one time is uniform over its area and thus the limitation of this second 
assumption does not apply. 

Any given part of a converging cylindrical or spherical shock front 
remains bounded by a certain wedge or cone, whose axis or vertex is 
respectively the centre of the cylindrical or spherical front. This axis or 
vertex is the origin at which the shock collapses. At a distance R from the 
origin the area of the wedge- or cone-shaped channel is proportional to R’ 
and hence the approximate relationship between shock strength and area 
becomes a (z, R)-relation. Near the origin it gives that the shock strength 
is proportional to R-/*, where K depends only on y. For the six cases 
studied by Butler the index —jK of this paper is found to be always 2°, 
of the index 2(a—1)/« of his exact solution. ‘This agreement is surprisingly 
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good. ‘The small discrepancy between the two solutions is due, as has 
already been pointed out, to the neglect of the reflected disturbances. By 
considering the etfect of these disturbances directly, it is found in §3 that 
the correction they cause is small due to an apparently fortuitous cancellation 
of terms. 

A check on the accuracy of this approximate theory away from the 
origin is provided by the work of Payne (1957). He uses a technique due 
to Lax (1954) to obtain, using an electronic computer, the flow behind 
converging cylindrical and spherical shock waves for various initial con- 
ditions. ‘The description of the motion of the shock is found to be in 
excellent agreement with the theory of this paper and is discussed in Payne’s 
publication. 

For problems involving sudden finite changes in area the work of 
Laporte (1954) applies. He determines the steady state solution for large 
area changes in a manner similar to that given in § 2 for small area changes. 


2. APPROXIMATE DESCRIPTION OF THE MOTION OF THE SHOCK 


If a shock moving uniformly along a channel encounters a small area 
change, the flow behind the shock and the shock itself are perturbed. As 
stated in the Introduction this problem has been solved by Chester. 
However the form of this solution at large time after the shock has passed 
the area change is easily obtained and is presented here. ‘The steady state 
configuration for the interaction of a shock and a small area change is shown 
in figure 1. 

Let p, p, u, a, M, with the appropriate suffix, be the pressure, density, 
fluid velocity, sound speed and Mach number in any of the six regions 
shown. ‘The incident shock strength z is defined as a pressure ratio, 


= (2.1) 
The flow variables behind the incident shock are given in terms of the 


shock strength and conditions ahead of the shock by the Rankine-Hugoniot 
equations which may be written in the form 


Pe (y—1)z+(y+1)’ (2.2) 
2p1 12 


2 12 


The fluid velocity u, is measured in the direction of the motion of the shock, 
and the velocity u, ahead of the shock is assumed zero. ‘The ratio of the 
specific heats of the ideal gas is y. 

Across the area change between regions 2 and 3, the continuity and 
Bernoulli equations together with the isentropic relation yield the first order 
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steady channel-flow equations, 
_yMz 6A 


a” 
Ps _ 5A 
Us 1 5A 
— =1+ (2.7) 
ly 
4 
! 
| 3 2 
5 
| 
2 
LE 
1 
6 
; 
! ' 
A At+SA DISTANCE 


Figure 1. A shock wave separating regions 1, 2 is incident on a small change in the 
area of achannel from A to 4+6A. The resulting transmitted shock separates 
regions 5,6. The fluid initially in the area change is set in motion by the 
shock and separates regions 4, 5. A small reflected disturbance separates 
regions 3, 4, and moves along the channel in the same direction as the shock, 
if the flow behind the shock is supersonic. The shape of the channel is 
shown at the bottom of the figure. 


The division between regions 4 and 5 separates the fluid into the parts 
initially on either side of the area change. After the disturbances have left 
this part of the fluid there will be no pressure or velocity differences between 


F.M, U 


| 
> 
| 
4 
f 
| : 
M 
4 


290 R. F. Chisnell 


these two regions, though it will be seen that a density change is left. Thus 


Pa=Ps = Us. (2.8) 
On passing the area change the shock encounters the same conditions 
as before but nevertheless its strength changes to z+6z due to the area 
change. The flow in region 5 is described by a set of equations similar to 
equations (2.1) to (2.4) with z replaced by z + 62 and the suffix 2 by suffix 5, 
it being remembered that conditions in regions 1 and 6 are identical. 
Combining all the above pressure and fluid velocity equations the 
following relations between the pressures and velocities in regions 3 and 4 
are obtained, again retaining only terms of order 5A, 


ps _ 1, 88, BA 


Ps 1A’ (2.9) 
My _ (y+1) 


But the pressure and velocity jumps across a sound wave are related 
directly, 
Ps 
In figure 1 this reflected disturbance is shown moving upstream with the 
shock, that is, the flow behind the shock is supersonic. The results of this 
section still hold if the flow is subsonic. 

Substitution of (2.9) and (2.10) into (2.11) leads to a relation between 
the first order quantities 54 and dz. This relation involves .i/, and M, 
which differ by a term of order 5A/A; thus, letting 54 > 0, both may 
be replaced by (2.4). There follows, after some reduction, 


1dA 1 1 (y +1) 


1/2 
‘Lage * 

An equivalent formula is given by Chester (1954) with the right-hand 
side of the equation written as {(z—1)K(z)}"1, where K(z) is expressed in 
terms of the Mach number of the fluid behind the shock, the Mach number 
of the shock referred to the fluid in front of the shock, and the relative 
Mach number of the shock to the fluid behind it. | Chester observes that 
K(z) is a monotonic decreasing function of shock strength with only a 
small total variation. For y = 1-4 this variation is from 0-5 for weak shocks 
to approximately 0-394 for strong shocks. 

It is to be emphasized that the derivation of (2.12) by matching steady 
state flows in regions 2 to 6 will give the solution only for large time. The 
work of Chester, however, is an exact solution to the problem depending 
only on the linearization with respect to the area change. Chester shows 
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that the shock strength averaged at one time over the shock front depends 
only on the area of the channel. Thus after the area change the average 
shock strength z+6z does not alter and is given by the above analysis. 
The only subsequent changes in the shock strength are local as the shock 
tends to a uniform plane shock. 

Integration of (2.12) will give a shock-strength/area relation applicable 
to channels of continuously varying area with finite total changes. When 
the area change is small this result gives the average strength of the shock 
after the area change, the shock then tending to become uniform with this 
average value as its strength. For large area changes however, the result 
gives only an approximation to the average shock strength. This is because 
disturbances reflected by the shock are not negligible when the area change 
is large. ‘These disturbances move into a non-uniform flow and thus 
generate further disturbances, some of which catch up the shock and alter 
its strength. It is the neglect of these disturbances that renders this 
shock-strength/area relation only an approximation to the average shock 
strength. The effect of these disturbances is considered in $3 for the 
cases of converging cylindrical and spherical shocks near the origin. 

The integration of (2.12) gives 


A f(z) = constant, (2.13) 
where 
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The function has been tabulated by Payne for y = 1-4 and is given in 
table 1. Given the value of the shock strength on encountering an area 
change in a channel and the ratio of the areas at the ends of the variable 
area section the strength of the shock on emergence from this section is , 
obtained by inverse interpolation in table 1. 

A much simpler, though not so accurate, relationship between shock 
strength and area is obtained by noting that A(z) has a small variation for 
quite large variations in shock strength. Thus replacing the right-hand 
side of (2.12) by {(z—1)K(sz)}" and treating K(z) as a constant, one finds 
that 


A(z—1)"* = constant. (2.15) 
U2 
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For y = 1-4 the values of K(z) are also given in table 1. It is suggested 
that an estimated mean of the values of K(z) corresponding to the shock 
strengths at the ends of the variable area section be used in (2.15). 

The application of the above description of a shock wave moving in a 
channel of varying area to converging cylindrical and spherical shocks is 
obtained by considering channels with areas proportional to R’, where 7 is 
1 and 2 respectively, and R is the distance of the shock from its axis or 
point of symmetry. 


F(z) K(z) z 10-* x f(z) K(z) 
1:00 0 0-500 00 jie 0-301 09 0-409 73 
1-05 0-019 119 0-493 37 14 0-453 05 0-408 04 
1-10 0-078 499 0-487 65 16 0-643 74 0-406 69 
1:15 0-181 03 0-482 65 18 0-876 10 0-405 57 
1-2 0-329 42 0-478 27 20 1:1529 0-404 65 
0-773 94 0-470 94 25 2-056 3 0-402 88 
1-4 1-431 0 0-465 06 30 3-291 7 0-401 62 
1°5 2-318 0 0-460 25 35 4-893 6 0-400 68 
1:6 3-450 6 0-456 26 40 6°894 1 0-399 95 
1:8 6:5109 0-449 99 45 9-323 1 0-399 36 
2:0 10-719 0-445 29 50 12-209 0-398 89 
2°5 26:867 0-437 40 60 19-457 0-398 15 
3-0 52-063 0-432 39 70 28-838 0-397 62 
3°5 87-418 0-428 81 80 40-539 0-397 21 
4-0 133-92 0-426 07 90 54-732 0-396 88 
4°5 192-48 0-423 85 100 71-582 0-396 62 
5 263-94 0-422 00 150 200-88 0-395 82 
6 448-72 0-419 03 200 417-44 0-395 41 
7 694-18 0-416 70 250 735-98 0-395 16 
8 1005-8 0-414 81 300 1169-6 0-394 99 
9 1388-6 0-413 24 350 1730-1 0-394 87 

10 1847-4 0-411 90 400 2428°5 0-394 78 


Table 1. The functions f(z) and K(z) relate to equations (2-14) and (2:15). y = 1-4. 


Thus for converging symmetrical shocks (2.13) becomes 
R f(z) = constant. (2.16) 
The simpler form (2.15) can be similarly applied. 

For weak shocks f(z) behaves like (s—1)? giving the familiar acoustic 
result that the strength of the disturbance varies inversely as the square 
root of the area of the disturbance. 

In the neighbourhood of the origin where is large f(z) behaves like z'", 
where K is the limiting value of K(z) for very strong shocks and is given by 


1 1 y 


Thus for cylindrical and spherical shocks near the origin the strength is 
proportional to R-*: and R-**: respectively. ‘Table 2 gives the exponents 
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—jK for various values of y, together with the corresponding exponents 
2(a—1)/x calculated by the exact theory of Butler. It is seen that the 
agreement is very good. For these symmetrical shocks it is to be remem- 
bered that the second limitation on the theory of this section, namely that 
it gives only an average shock strength, does not apply. The only source 
of error is the neglect of the effect on the flow behind the shock of the 
reflected disturbances. This will now be considered. 


Cylindrical Shock (j = 1) Spherical Shock (j = 2) 
—K 2(a—1)/a —2K 2(a—1)/x 
this paper Butler this paper Butler 
y= — 0-326 223 —0-322 441 — 0-652 447 —0-641 513 
—0-394 141 — 0-394 589 —0-788 283 — 0-788 728 
y=3 —0-450 850 — 0-452 108 —0-901 699 —0-905 385 


Table 2. The values of —K or —2K and 2(a—1)/a are the distance exponents 
describing the behaviour of the shock strength, for cylindrical and spherical 
shocks near the origin, on the theory of this paper and the exact solution of 
Butler respectively. 


3. DISCREPANCY BETWEEN THE PRESENT THEORY AND THE SIMILARITY SOLUTION 


It was seen in the last section that the description of symmetrical shocks, 
obtained from Chester’s channel result, is in close agreement with that 
given by the similarity solutions of Guderley and Butler. Such discrepancy 
as there is may be indirectly attributed to the disturbances generated by 
the motion of the incident shock. As these disturbances move through the 
non-uniform medium behind the shock, their strengths change and further 
disturbances are reflected which move through the fluid in the same 
direction as the shock. The strengths of these re-reflected disturbances 
are examined in this section. In particular the modification to the motion 
of the shock due to some of the re-reflected disturbances merging with it 
is discussed. 

A complete account of the effect on the shock of disturbances merging 
with it from behind is not given, there being three limitations to this work. 
The first is that the effect on the shock of higher order disturbances—such 
as the fourth order disturbance generated by the third order or re-re-reflected 
disturbances—is neglected. ‘The second limitation to the work of this 
section is an assumption concerning the description of the flow behind the 
shock. It is known that the shock strength is singular at the origin, giving 
the Mach number of the flow behind the shock the limiting value 


M= (3.1) 


and the density ratio across the shock the limiting value (y+ 1)/(y—1). 
The reflected disturbances and their associated re-reflected disturbances 
which merge with the shock are assumed to travel in a flow having this 
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limiting Mach number and density. Due to the violent increase in shock 
velocity as it approaches the origin it is expected that the major part of the 
disturbance which catches up the shock before it reaches the origin is 
generated in the neighbourhood of the origin. However there will be 
some disturbance generated away from the origin which in the early part 
of its motion towards the shock moves in a region where the above limiting 
flow values are not applicable. ‘The third assumption is that elements of 
the reflected and re-reflected sound waves travel without changing their 
strength. ‘The author has attempted to improve the description of the 
re-reflected wave given in this section by considering the various causes 
of changing strength of the elements of the sound waves, but the attempt 
has not been successful. 

There are three types of re-reflected disturbance, which together 
produce only a small effect on the shock. One of these does itself produce 
only a small effect, but each of the other two has itself a considerable effect 
on the shock; it is the nearly complete cancellation of these three types of 
re-reflected disturbance that renders the description of the shock in §2 so 
close to the exact similarity solution in the neighbourhood of the origin. 
The completeness of this cancellation can be illustrated by reference to 
table 2, where it is seen that, in the case y = 3 for example, the exponent 
describing the motion of the shock deduced in § 2 agrees with the similarity 
solution to two significant figures. If, however, just one of the two more 
significant types of re-reflected disturbance is considered the exponents 
differ by one unit in the first significant figure. 

The reflected disturbances are first considered. It was seen in §2 
that there are two types of disturbance generated by the shock, a sound 
wave and an entropy variation. ‘The disturbances produced by the shock 
as it moves through an elementary change in the channel area are shown in 
figure 1. An element of the sound wave separates regions 3 and 4, and a 
small entropy variation exists between the neighbouring fluid particle 
paths separating regions 4 and 5. ‘The entropy variation is represented 
as a change in density between these regions, which have the same pressure, 
and will be referred to as a weak contact discontinuity. ‘he strength of 
these disturbances on reflection may be obtained from the analysis of the 
last section. ‘The strength of an element of the sound wave, defined as a 
pressure ratio, is given by (2.9) as 
* A 
after (3.1) has been used to eliminate M,. The strength of the weak contact 
discontinuity separating regions 4 and 5 in figure 1 may be found by the 
following simple argument. When a ‘strong’ shock moves through an 
area change 6A there will be no change in the density ratio across the shock, 
which has the limiting value (y + 1)/(y—1), although there will be a change 
in the pressure ratio across the shock (dz/dA)5A.__In the absence of entropy 
changes the density ratio would be correspondingly changed by a factor 
1+ 1/y3(dz/dA)bA. It is this change in density ratio that must be balanced 
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by the reflected weak contact discontinuity. Hence 


1 dz : 
waa 5A. (3.3) 

The coefficient dz/dA is given by (2.12), which for strong shocks assumes 
the form 

1 dA 1 
(3.4) 

The generation of an element of the re-reflected wave is due partly to 
the interaction of elements of the two forms of reflected disturbance and 
partly to the motion of these elements through changes in area. This 
interaction is shown in figure 2 to be taking place at a point in the channel 
where the area is A,,. An element of the sound wave and a weak contact 
discontinuity which arrive at A,, simultaneously are shown to have been 
generated by the shock at points of the channel where the areas are A, 
and A,, respectively. 


TIME 


Ac A, A. Ay DISTANCE 


Figure 2. An element of the sound wave reflected by a shock wave is generated at a 
channel area of Ag. It interacts with a weak contact discontinuity at a channel 
area A,,, the contact discontinuity having been left behind in the fluid by the 
shock at a channel area A;.. The re-reflected disturbance resulting from this 
interaction and the passage of the two disturbances through a small area 
change at A,,, merges with the shock wave at a channel area A yy. 


e 
IS 
‘t 
f 
r 
) 
| 
y 
) 
| / 
| ' 4 
1 
/ 
| i 1 
1 
| 


296 R. F. Chisnell 


The strength of the element of the sound wave on arrival at A, is assumed 
to be unaltered from its value upon being generated at A,. This is given 
by (3.2), which may be rewritten, using (3.4), as 
dA, 


2y 
where = { - Kt. 3.6 
The strength of the elementary contact discontinuity will not be altered 
by the constant entropy interactions occurring behind the shock. Hence 
its strength on arrival at A,, is the same as its strength on generation at A,,, 


which is given by (3.3) as 


l+y, (3.5) 


(3.7) 


again using (3.4). 

It is not intended to present here the analysis of the interaction which 
leads to the generation of a typical element of the re-reflected wave. The 
method is that used in §2, in which steady state solutions are matched 
across the various disturbances present. A very similar interaction is 
studied in detail in a previous publication by the author (Chisnell 1955). 
The strength of the re-reflected element in excess of 1 is found to be 


V1 K dA, K dA, dA, 


(2-yP(y+ 
BY? 7 +2)— Di 

The three terms in (3.8) represent respectively re-reflections due to 

(i) the interaction of an element of the sound wave and a contact discon- 

tinuity, (ii) the motion of the contact discontinuity through anareachangedA,, 

and (iii) the motion of the element of the sound wave through this area 

change. ‘To compare the terms in (3.8) it is necessary to express 6A, 5A,,, 

5A, 6A,,, 5A, 5A, in common form. To do this, new variables A,,, the 

area of the channel where the re-reflected disturbance merges with the shock, 
and A given by 


(3.8) 


(3.9) 


Ay\" 1 K 1 


are introduced. ‘Iwo approximate relations between A, A, A, and A,, 
are formulated, only two of these variables being independent. ‘These two 
relations express firstly that the element of the sound wave and the weak 
contact discontinuity arrive at A,, simultaneously and secondly that the 
re-reflected disturbance merges with the shock at A,,. The relations 
are obtained by integrating the inverse of the speeds of the various 
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disturbances with respect to distance, hence obtaining simple equations 
expressing the information that given disturbances arrive at specified 
points simultaneously. Behind a strong shock the fluid velocity and the 
sound speed are proportional to the shock velocity, the constants of 
proportionality being 2/(y+1) and »/{2y(y—1)}/(y+1) respectively. All 
these velocities are singular, being like z¥? or A-*” from (3.4). The distance 
along the channel is proportional to A", 7 being 1 or 2 for the cylindrical 
or spherical shock. Hence the integral of the inverse of speed of disturbance 
with respect to a distance will be proportional to A”. The two required 
relations are therefore- 


y+1 


ytl n ytl n_ qn 
(AX A’) V {2y(y - 1)} (As Aj) + 2+ V {2y(y 1)! (Aj; A‘). 


(3.13) 
In terms of the variables A,,, A, (3.8) becomes 


By forming an integral from (3.14) with respect to A, one finds that the 
total strength of the re-reflected disturbances which merge with the shock 
at the area A,, is given by 


_, 


1 (3.16) 


Ky | (1-ya)ys top (1—y,). (3.17) 


og 
4yn (y¥s— 74) yn (y¥s— 4) n 
The integrations with respect to A are between the limits 0 and y,, repre- 
senting re-reflected disturbances formed a long way behind the shock and 
just behind it, see (3.11). 

A description of the motion of the shock allowing for the effect of the 
re-reflected disturbances merging with it may be obtained by modifying 
the work of §2. In figure 1, region 2 will now be split into two parts by the 
re-reflected disturbance. Combining (2.16) and (2.17) the original shock- 
strength distance relation may be written as 


sR* = constant. (3.18) 


Allowing for re-refiections this becomes 


zR#*® — constant, (3.19) 


where 1 y—-1\) (3.20) 
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Cylindrical Shock (j = 1) Spherical Shock (j = 2) 

(i) (ii) (iii) (i) (ii) (iii) 
| +0-0072 | —0-049 | +0-005 2 | + 0-013 —0:085 | +0-0075 
y= % | —0-00019} -+0-058 | —0-00052} —0-000 33} +0-099 | —0-000 26 
y= | —0-0017 | +0-095 | —0-0013 | —0-0029 | +0-16 —0-002 0 


Table 3. A direct consideration of the re-reflected wave gives the estimates (i) for 
n, which are compared with (iii) the values of » required to bring the theory 
of this paper into exact agreement with the similarity solution. The entries 
(ii) give the largest single contribution to the values of 7 given in (i). 


The cancelling between the three types of re-reflected disturbance is 
shown in table 3, which gives (i) the values of 7 obtained from (3.17), and 
(ii) the largest of the three terms in this expression for 7. For the various 
cases considered it is seen that 7 is only between } and 54 of the largest 
term. Also given in this table are (iii) the values of » required to give exact 
agreement with the similarity solution. It is seen that the values (i) are 
of the right order of magnitude to explain the small discrepancy between 
the work of §2 and the similarity solution. Thus the cancellation between 
the three types of re-reflected disturbance explains the very close agreement 
of the work of this paper with the similarity solution. The cancellation of 
the terms does however limit the accuracy with which the discrepancy 
between the two theories may be directly estimated, consistent with the 
assumptions of this section. As the third assumption, that the elements 
of the sound waves move with unchanged speed, applies only to the first 
and third term of (3.14), it is probably the most critical of the assumptions. 
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Diffusion from a fixed source at a height of a few 
hundred feet in the atmosphere 


By J. S. HAY and F. PASQUILL 
Meteorology Section, Chemical Defence Experimental Establishment, Porton, Wiltshire 


(Received 19 January 1957) 


SUMMARY 


The results of measurements of the vertical distribution of 
airborne particles, released usually at a height of 500 ft., and 
sampled for periods of about 30 minutes at downwind distances 
of 100, 300 and 500 metres, are presented and discussed. 

At all distances the frequency distributions of the particle 
elevation with respect to the site of release are closely similar in 
shape and size to the frequency distribution of wind inclination at 
the site of release. ‘This is interpreted as showing that, despite the 
uncorrelated effects of small eddies, high correlation was maintained 
in the motion of a particle, for periods of 1 minute or more, by 
the dominant action of persistent eddies which contributed heavily 
to the turbulent energy. In contrast to this result, the wind 
observations showed that the Eulerian auto-correlation coefficient 
(R,) fell to about 0-2 in 10 seconds. 


1. INTRODUCTION 
The turbulence which occurs in the atmosphere clear of the ground 
is often approximately steady and homogeneous. Under these conditions, 
the diffusion of particles released at a fixed point is specified by the well- 
known ‘Taylor (1921) relation 


= | Re (1) 


Here z is the vertical displacement of particles from their mean vertical 
position after travelling a time 7, w’ is the deviation from the mean of the 
vertical component of the velocity of the particles, and R, is the correlation 
coefficient between the vertical velocity at one instant and the velocity of 
the same particle at a time € later. Rg is usually referred to as the Lagrangian 
correlation coefficient. In practice Rz is very difficult, if not impossible, 
to measure effectively and in problems of atmospheric diffusion progress 
has depended (see Sutton 1932, 1934) on specifying intuitively or empirically 
some general relation between R;z and measurable quantities. One empirical 
approach, however, which does not seem to have been exploited to any 
extent, is to use experimental measurements of diffusion to obtain some 
indications of the properties of the Langrangian correlation coefficient in 
comparison with those of the corresponding, measurable, Eulerian coefficient 
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Ry, that is, the correlation between the velocity of the air at a point at one 
instant and that at the same point at time ?¢’ later. The simplest possible 
type of comparison would be to evaluate 


1 
| dt'dt 
and compare it, as a function of 7, with the corresponding double integral 
from (1). After commencing the present investigation the writers’ attention 
was drawn to some recent work by Mickelsen (1955) in which precisely this 
line of approach has been followed in studying diffusion in homogeneous 
turbulence in a wind tunnel. 

The present paper describes experiments in which the vertical distri- 
bution of airborne particles, released in the atmosphere at a position remote 
from the ground, was examined at various distances downwind up to 
500 m, using simple lightweight apparatus carried on the cables of captive 
balloons. Simultaneously, measurements of the fluctuations in the 
inclination of the wind at the site of release were made, in order effectively 
to specify w’ and R,, using a balloon-borne gustiness recorder described 
previously by Jones (1956). 


2. EQUIPMENT AND TECHNIQUE 

The particles used in these experiments, spores of Lycopodium, were 
chosen because they are large enough to impact with reasonable efficiency 
on cylinders of convenient diameter, and can thus be collected from the air 
without the aid of elaborate sampling apparatus. This method of sampling 
was used successfully by Gregory (1951) in wind tunnel experiments. 
As these spores come from a type of moss usually found in mountainous 
districts, it was unlikely that they would be present naturally in the 
atmosphere over downland near Salisbury, a fact confirmed during the 
course of the experiments. Their characteristic size (304 diameter), 
shape and colour enable them to be distinguished without undue difficulty 
from other particles suspended in the atmosphere, while their terminal 
velocity (about 2cm/sec) seems low enough not to affect their scattering in 
a turbulent atmosphere. 

A dispenser, shown in figure 1 (plate 1), emitted the particles in a 
quasi-continuous manner, at an average rate of 3gm/min. The spores are 
gravity fed from a hopper (of dimensions 17 x 6x1 in.) onto a toothed 
wheel rotating in a Perspex block, from which they emerge into the air 
and onto the blades of a small high-speed fan which serves to break up 
aggregates of spores. ‘The spores are emitted in discrete puffs but the 
interval between puffs is only about } sec. In practice, the dispenser was 
mounted on the cable of a large captive balloon, usually at about 500 ft. 
above the ground, and operated for 30 min or so at a time. Power for the 
motor driving the toothed wheel and fan was supplied from a source on 
the ground via lightweight electric cables. 
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The apparatus for collecting the particles consisted simply of 5 cm 
lengths of } in. glass tubing, of which some 4 cm were covered with an 
adhesive similar to that described by Gregory (1951). Cylinders of this 
size are easily handled and examined and, according to Gregory, their 
collection efficiency* for Lycopodium spores in winds of more than 3-3 m/sec 
exceeds 40% and does not vary appreciably with wind speed. Throughout 
each emission, the adhesive cylinders were held into wind by small vanes 
spaced along the cable of a small captive balloon so that the central vane 
was about level with the dispenser. This balloon was positioned downwind 
of the balloon carrying the dispenser, at a distance of 100, 300 or 500m, 
and was moved as necessary to allow for any substantial changes in wind 
direction during the emission. 

In the subsequent assessment of the samples collected, twenty sections 
of each cylinder were examined by a microscope having a field of view about 
():12 cm square and the number of particles in each section was counted. 
These numbers were found to follow a Poisson distribution so that the 
mean value, which was taken to be a measure of the deposit on a cylinder, 
had a standard error which was easily calculated. 

Mounted immediately above the dispenser (see figure 1) was the 
instrument for the simultaneous continuous recording of the fluctuations 
of wind inclination and of mean wind speed (Jones 1956). Thermistors 
were also mounted on the cable and connected to recording apparatus on 
the ground in order to give the temperature difference between the level 
of the dispenser and a point 200 ft. below. The temperature difference 
between 23 ft. and 4 ft. above the ground was obtained from permanent 
apparatus located nearby. To allow corrections for the vertical movement 
of both balloons to be applied, observations of the elevations of the dispenser 
and central sampling vane were made by means of a theodolite. The time 
interval between observations was usually 10 sec in the case of the dispenser 
and 1 min in the case of the central sampling vane. 


3. REDUCTION AND PRESENTATION OF OBSERVATIONS 
(a) Particle distributions 


Of the 19 experiments which were carried out, it was evident that in 
9 cases the cloud of particles had not been adequately sampled, due to 
insufficient spacing or incorrect downwind positioning of the sampling 
cylinders, and these were not considered further. In the remaining 
10 cases, to facilitate comparison with the measured wind inclinations 
at the site of release, the vertical distribution of particles was expressed 
in terms of their angular elevations with respect to the dispenser. The 
particle deposit at each level was taken to be representative of the interval 


* The collection efficiency is the number of particles impacting on the cylinder, 
expressed as a percentage of the number which would have passed through the same 
cross-sectional area if the cylinder had not been there. For a given particle size, 
the collection efficiency depends on the diameter of the cylinder, the speed of the 
air stream and the type of adhesive used. 


=¥ 

% 

at 

Fi 

: 

2 


302 Jj. S. Hay and F. Pasquill 


centred on the cylinder and extending mid-way to the adjacent sampling 
levels. Cumulative deposits, obtained by addition of successive individual 
values, and expressed as a percentage of the total, are plotted against 
elevation on probability graph paper in figures 2 and 3. 
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Figure 2. Distributions of particle elevation (x) and wind inclination (©) for 
experiments 4, 7, 8, 16 and 17. Ordinate scale is given for experiment 4. 
Scales in other cases are similar with zero lines indicated by arrows. 


The observed spread of particles contains a contribution arising from 
the relative vertical motion of the dispenser and sampling array due to 
balloon movements. This contribution depends on the difference, AH, 
between the height of a datum point on the sampling array at ‘time 7 and 
the height of the dispenser at time T—t, where t¢ is the time of travel of 
the particle. It seems unlikely that there could have been a high degree 
of correlation between AH and the true vertical displacement of the particle 
with respect to its point of origin. This supposition was confirmed by 
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examining the relation between AH and ¢, where ¢ is the wind inclination 
at the dispenser position at time 7—t, for experiment no. 8, which 
represents the worst case in the sense that the variations of AH relative 
to the apparent vertical displacements of particles were largest. The 
usual statistical rule for deriving the variance of the difference of two 
independent quantities may therefore be applied. If o? and oc? are 
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Figure 3. Distributions of particle elevation (x) and wind inclination (©)/ ‘for 
experiments 12, 13, 14,15 and 19. Ordinate scale is given for experiment 12. 
Scales in other cases are similar with zero lines indicated by arrows. 


respectively the variances of AH and the apparent vertical displacements 
of the particles, then O°, the variance of the true vertical displacements, 
is gi 
is given by = 
In practice, these quantities were expressed as elevations with respect to 
the dispenser. The values of AH were obtained from the theodolite 
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observations, the time ¢ being estimated from the wind speed measured 
over periods of a few minutes. 


‘b) Distributions of wind inclination 

From the gustiness records of these same 10 experiments, values of the 
wind inclination ¢ were read off to the nearest }° at 2}-sec intervals. Because 
of the cramped nature of the records, it was frequently possible to estimate 
only upper and lower limits for the value of ¢, and in such cases the mean 
of these limits was taken as the appropriate value. Corrections for 
deflections of the wind-inclination vane, which arose from the movement 
of the supporting balloon rather than from the turbulent motion of the air, 
were then applied as follows. 

If the elevation « of the wind-inclination vane as observed by the 
theodolite is small, and changes by an amount Aq in a time A¢ (usually 
10 sec), the spurious component of wind inclination over this time interval 
is approximately ¢, = —/Aa/uAt where / is the horizontal distance from 
the theodolite to the vane (about 3000 m) and u is the indicated horizontal 
wind speed over the interval. The error is principally due to the neglect 
of the horizontal components of the eddy motion and the motion of the 
balloon. In evaluating the expression, it was assumed that values of u 
averaged over a period of 2 or 3 min applied to the component intervals At. 
Instantaneous values of ¢, at times corresponding to the readings of ¢ 
were then estimated to the nearest $° from these short period means by 
graphical interpolation, and used to correct the values of ¢. 

The corrected values of wind inclination, that is, 6—¢,, from each 
experiment were formed into a frequency distribution based on class 
intervals which were one degree wide and centred on integral values of 
inclination, the frequencies of the intermediate half-degree values being 
divided equally between adjacent intervals. For comparison with the 
particle distributions, cumulative frequency distributions were formed for 
the range of elevations covered by the sampling cylinders. This procedure 
necessarily excluded the occasional extreme values of inclination, which 
were mainly positive in sign, that is, indicative of upward moving cu*rents. 
These distributions are shown together with the corresponding particle 
distributions in figures 2 and 3, the cumulative frequencies being plotted 
against the upper limits of the respective class intervals. 


(c) Meteorological and other conditions 

The dates and times of the various experiments, and the distances and 
heights involved, are given in table 1 together with wind speed and local 
and surface temperature gradient measurements. Mean wind speeds 
at the height of release ranged from 6-1 to 9-5 m/sec and the 
temperature gradient in the lower layers ranged from slight lapse to 
slight inversion. 
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at a height of a few hundred feet in the atmosphere, Plate |. 
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Figure 1. The particle dispenser and wind inclination instrument on balloon cable. 
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Experiment Number 4 7 8 12 13 
Date (1955) 16 Sept. 26 Sept. 26 Sept. 4+ Oct. + Oct. 
Time (G.M.T.) 0828-0844 | 1436-1506 | 1606-1636 | 1521-1601 ; 1636-1716 
Distance of sampling 

m 100 100 100 300 300 
Height of dispenser 

above ground (m) 143 147 255 140 143 
Wind speed (m/sec) 6°7 8-0 7:4 7:7 7:33 
Local temperature 

gradient 0-90r 0-661 0-621 
Surface temperature 

gradient (23’—4’) 

CF) 0-2 0 -Q-4 0-6 
Cloud amount and As Se 5000 Sc Sc 5000 | Se 5000 

height (ft.) # Se 2500 % Cu 3000 3500 = Cu 3000 | 2 Cu 3000 
Experiment Number 14+ 15 16 17 19 
Date (1955) 13 Oct 14 Oct. 17 Oct. 17 Oct. 25 Oct 
Time (G.M.T.) 1522-1552 | 1425-1455 | 1151-1221 | 1430-1500 | 1517-1557 
Distance of sampling 

(m) 300 300 500 500 300 
Height of dispenser 

above ground (m) 145 154 150 151 129 
Wind speed (m sec) 6°3 6-1 3 6°8 9°5 
Local temperature 

gradient 0-82] 0-861 0-870 0-761 
Surface temperature 

gradient (23’-4’) 

Cloud amount and 2 CiCs Se 2500 AcAs CuSe Sc 

height (ft.) 13 Cu 3—-% CuSc 3000 3000 

2500 3000 


Table 1. Summary of meteorological and other conditions. 


The local temperature gradient 


was obtained from measurements over a height interval of 200 ft. near the dispenser. 
I’ = dry adiabatic lapse rate. 


4. ‘THE OBSERVED RELATION BETWEEN THE DISTRIBUTIONS OF PARTICLE 


ELEVATION AND WIND INCLINATION 


The frequency distributions of particle elevation and wind inclination 
collected in figures 2 and 3 show three features of interest. Firstly there is 
the relation between the median values of the two distributions. For seven 
of the experiments there is agreement to within 1°, in one experiment the 
particle elevation is higher by 3°, while in the remaining two, it is lower 
by 2°. The mean of the differences in median values, in the sense wind 
inclination minus particle elevation, is +0-2° with a standard deviation 
of 15°. In other words, there is no significant indication of a declination 
of the trajectories of the particles with respect to the initial air flow. However, 
it is noteworthy that this value +0-2° is very close to that which would 
follow from the expected terminal velocity of 2 cm/sec in a wind of about 
7 m/sec. In individual experiments the gravitational fall is evidently 
negligible in comparison with the unknown effects of incomplete sampling 
(in a statistical sense), and of systematic features such as topographical 
control on the general flow pattern or variations of eddy structure with 
height. 

F.M. 
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A second feature is that the graphs in figures 2 and 3 are roughly linear, 
indicating an approximation to Gaussian form in the distributions of both 
wind inclination and particle elevation. ‘This is not a particularly novel 
conclusion concerning the wind data—except in so far as it refers to wind 
inclinations at an unusually elevated site—for it is widely thought that 
turbulent components will usually be so distributed, apart from anisotropic 
motions associated with well-developed convective action. In this latter 
respect it is to be borne in mind that the restriction of our comparison of 
wind and particle distributions to the height range covered by the sampling 
apparatus has necessarily excluded the occasional extreme inclinations, 
most of which are a reflection of this anisotropic feature. The main interest 
is rather in the approximately Gaussian form of the particle distributions, 
for here it is to be remembered that the measurements are not made at a 
single level of say 500 ft., but over a layer which may extend between 200 
and 800 ft. and which will contain systematic variations of wind and eddy 
structure with height. Inspection of the individual distributions does not 
suggest any consistent departure from Gaussian form and it is interesting 
to note that there is no outstanding difference in this respect between results 
obtained at different distances from the source. 

The third and most interesting feature is the similarity in the dimensions 
of the wind and particle distributions, though before any realistic comparison 
can be made it is necessary to correct for the relative vertical displacements 
of the dispenser and the sampling array, as discussed previously. Values 
of o,, the standard deviation of the apparent elevations of the particles, 
were extracted from the distributions in figures 2 and 3 by reading off the 
elevations at 2 and 98%, that is, those enclosing 4:12 times the standard 
deviation in a Gaussian distribution. A similar process was used to find 
o,,, the standard deviations of wind inclination. ‘The standard deviations, 
o,, of the elevations due to the relative vertical displacements were computed 
directly. ‘The values of o, and o,, the resulting values of o,, and the 
corresponding values of o,,, are collected in table 2 below. Examination 
of the statistics of sampling indicates that the standard errors in the individual 
values of particle deposit range from a few per cent for high deposits to 
greater than 20%, for the very low deposits. The net error in o,, due to these 
errors is estimated to be of the order of 3%. In the case of o,, the main 
error is due to the uncertainty in the correction applied for the movement 
of the balloon. An upper limit to this error has been estimated for one case 
(no. 17) and found to be of the order of 4%. 

With the exception of nos. 4 and 19 there is a remarkable degree of 
similarity between the values of o, and o,,, irrespective of the distance of 
travel and the conditions of wind and atmospheric stability prevailing. 
The odd result in the case of no. 19 is clearly associated with certain subsidiary 
peaks in the particle distribution, which occur at a height well below the main 
peak and for which there appear to be insufficient corresponding large 
negative values of wind inclination. No. 4 was carried out in pre-frontal 
conditions with drizzle setting in, and is a special case in that an abrupt 
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reduction in gustiness occurred about half-way through the experiment. 
The observation of a distribution of particles which is narrow relative to 
the distribution of wind inclination could be due to unrepresentative 
sampling, since the positioning of the sampling array with respect to the 
cloud of particles was probably controlled more successfully in the second 
half of the experiment than in the first. ‘This idea is supported by the very 
high particle deposits that were observed over a sampling time of only 
16 minutes. 


Dist. Particle elevation Wind re 
AT°F 
Expt. of inclin- u 
(23 
no. travel Ca o, | ationo, | (m/sec) 4 ft.) 
(m) (degrees) (degrees) 
4 100 1-2 2-1 1-7 0°55 6-7 
7 100 1:3 3-9 3:7 3-0 1:23 8-0 0-2 
8 100 2-2 4-1 3°5 2°8 1-25 7:4 0 
12 300 0:7 5-0 49 4-6 1-07 77 —0-4 
13 300 3-9 3-9 3°8 1-03 7:3 
14 300 1:0 4-2 4-1 4-0 1-03 3 —0-4 
15 300 0-2 4:9 4-9 4:9 1-00 671 —0-6 
16 500 0-5 3-0 3-0 3-2 0-94 7:3 —0:7 
17 500 0-3 2:9 2:9 2:7 1-07 
19 300 0:9 6-1 6:0 4-0 1:50 9°5 +0-2 


Table 2. Comparison of standard deviations of particle elevation and w.nd 
inclination. 


The original data have been independently examined by a statistician. 
That the distribution was approximately Gaussian, in one case which 
appeared doubtful, was confirmed by the application of a formal test. 
In the subsequent treatment, therefore, the data were assumed to be of 
Gaussian form. Variances of apparent particle elevation (o?) and of wind 
inclination (o2) were calculated directly from the observations. The 
values of o2 were then corrected as before to give o?. (The values of o, 
and o,, obtained in this fashion differed only slightly from those presented 
in table 2.) Owing to the large number of degrees of freedom involved, 
the equality of o? and o; could not conveniently be tested from the available 
statistical tables of variance ratio. Instead, a form of t-test was applied 
to the differences of the variances. Of the ten experiments, six (nos. 12, 
13, 14, 15, 16, 17) showed no significant difference, while for the remaining 
four (nos. 4, 7, 8, 19) the differences were found to be highly significant. 
However, only in no. 4 is o,, substantially less than o,,, and this is the 


experiment which has already been discussed as possibly anomalous. 


5. INTERPRETATION OF THE RESULTS IN RELATION TO EDDY STRUCTURE 


As stated in the Introduction to this paper, the fundamental object of 
the present investigation was to obtain indirectly some insight into the 
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Lagrangian features of the eddy structure, that is, variations which are 
related to time in the sense of following the motion, in comparison with 
Eulerian features, which are related to time at fixed points in space. From 
Taylor’s expression for the scattering of particles (equation (1)), it immediately 
and necessarily follows that a linear variation of standard deviation with 
time (in our notation, o, = o,,) is held as long as the Lagrangian correlation 
coefficient remains close to unity. For such a time, the particles travel in 
straight lines with the velocity imparted to them at the moment of release. 
It should be noted, however, that the departure from the linear variation 
of o,, is not sensitively dependent on the initial fall in the correlation 
coefficient. For example, if this fall is at an exponential rate, the value 
of a,,/o,, will still be as much as 0-86 when the correlation coefficient haz 
fallen to about 0-4. 

The present results, in showing equality between o,, and a, for distances 
up to 500 m, imply that the diffusive spread of the particles may be described 
by equation (1) with R; maintained at a relatively high value for at least 
1 min. This period of maintenance of Lagrangian correlation is considerably 
longer than that which is evident in the Eulerian or auto-correlation 
coefficient computed from corrected values of the wind inclination* at 
2}-sec intervals, as can be seen from the two cases quoted below in table 3. 
These figures are the means of three values obtained from analysis of three 
separate 8-min periods contained in each of the two experiments. 


| 


Time lag (sec) 24 10 15 20 30 


Correlation ( Expt. no. 13 0-69 0-40 0-19 0-09 —0-01 0-01 0-02 


coefficient Expt.n0.17 | 0-61 0-40 0:29 0:23 0-20 0-14 0-04 


Table 3. Values of auto-correlation coefficient for wind inclination. 


The simplest interpretation of the present measurements of Lagrangian 
correlation is that, despite the rapid Eulerian variations, the particles travel 
in virtually straight lines for substantial periods. It is clear from practical 
experience, however, that this interpretation cannot be literally correct. 
Indeed, it is obvious from visual observations of a puff of smoke moving 
with the wind that the trajectories of individual smoke particles are far 
from linear, and are in fact rapidly distorted mainly by turbulent motions 
on a scale smaller than that of the puff itself. 


* For the purpose of the arguments which are ultimately developed in the present 
paper it is probably unnecessary to give detailed consideration to the differences 
between the auto-correlation for wind inclination (R;(¢)) and that for the vertical 
component of eddy velocity (R,(z)). However, an analysis of this feature has 
confirmed that small values of R,(é), of the order of 0-1, are necessarily associat: d 
with similarly small values of R,-(z). 
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‘The explanation of the above paradox clearly lies in the form of the 
spectrum of eddy motion existing in the atmosphere. ‘The trajectory of 
a single particle may be regarded as composed of a series of harmonic 
oscillations. If the oscillations of low and high frequencies have respectively 
large and small amplitudes, the particle will tend to follow a path which 
contains small-scale irregularities about a large-scale relatively smooth 
curve. ‘The latter will be virtually straight over appreciable distances, 
and will therefore lead to high values of the Lagrangian correlation coefficient 
for the particle. 

A more general analytical discussion of this effect of the Lagrangian 
spectrum of turbulence has been given by Batchelor (1949), who demon- 
strated that the statistical dispersion of a single fluid particle is determined 
initially by all frequencies (of the particle velocity) in proportion to their 
energies, but that at greater distances from the point of release the effects 
of the lower frequencies become relatively more important. The point 
which is brought out by the present experimental results in the atmosphere 
is that, despite the uncorrelated effects of small eddies, high correlation 
may be maintained in the motion of a particle by the action of persistent 
eddies which contribute heavily to the total energy of turbulence. It may 
at first sight seem surprising that, in contrast to this result, Edinger’s (1951) 
measurements, with soap bubbles in an isothermal layer at a height of 
1000 ft., indicate a Lagrangian correlation coefficient falling to 0-5 in about 
7 sec. However, this is presumably due to the fact that Edinger’s observa- 
tions were made over a total time of only a few seconds, and therefore could 
not include the effects of any large persistent eddies, even if these existed 
in the conditions of his experiments. 

The difference in the persistence of correlation in the Eulerian and 
Lagrangian systems is presumably associated mainly with the fact that 
in the Eulerian measurements the pattern of eddy motion is being swept 
more or less rapidly past the fixed observing point. ‘The effect has previously 
been observed in wind tunnel turbulence, and, indeed, in an investigation 
similar in some respects to the present one, Mickelsen (1955) has recently 
compared the diffusive spread of a gas emitted from a point in wind tunnel 
flow with that calculated from equation (1) with the Lagrangian correlation 
coefficient replaced by measured Eulerian values. Equality in the observed 
and computed forms of o?/2v?T?, where v is the cross-stream component 
of eddy velocity, was found to occur when 7 (Lagrangian) was approximately 
3 times 7' (Eulerian). ‘This relation held for a range of values of o?/v?T? 
from 1 to 0-4. The data did not permit evaluation at lower values of 
o?/y®T? but the trend appeared to be in the direction of a rapid increase 
in the ratio of the ‘ Lagrangian’ to the ‘ Eulerian’ time. 

In contrast to Mickelsen’s experiments, ours do not extend into the 
region of non-linear relation between o? and #*, and a similar evaluation 
of the ratio of the ‘ Langrangian’ to the ‘Eulerian’ time for equality of 
the above forms is accordingly not possible. However, substitution in 
equation (1) of the auto-correlation data given in table 3 leads to a 
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(c, t)-relation which begins to depart significantly from linearity (i.e. 
o?/w”t? = ()-8) for ¢ between 5 and 10 sec, whereas in the observed values 
a similar order of departure from linearity is not evident even at 70 sec. 
In other words, this suggests that for equal values of o?/w’*t?, the ratio of 
the ‘ Lagrangian’ to the ‘ Eulerian’ time is not less than about 10:1. 

In the observational study of atmospheric turbulence there have recently 
been some quite novel measurements made by Gifford (1955), in an attempt 
to compare the Lagrangian and Eulerian characteristics directly at a height 
of 300 ft. Gifford uses the trajectories of ‘no-lift’ balloons from which to 
extract vertical velocities of a Lagrangian nature, and compares these on 
an energy spectrum basis with Eulerian values from a fixed installation. 
The interesting feature of these results is the suggestion that the Lagrangian 
spectrum is generally displaced to the /ow frequency side of the Eulerian 
spectrum, with the respective frequencies of peak energy in the ratio of 1:3 
approximately. 

It seems clear, therefore, that our conclusion regarding the difference 
in Lagrangian and Eulerian persistence of eddy velocity is not inconsistent 
qualitatively with the results of Mickelsen and Gifford. 
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Corrigenda to Drift” 


By M. J. LIGHTHILL 
Department of Mathematics, University of Manchester 


(Received 15 March 1957) 


In the author’s paper ‘“ Drift”? (Lighthill 1956) there are mistakes which 
need correction. Equation (28) should read 


m 


since the definition of the disturbance potential ¢ (p.35) states that the 
full velocity potential is U(v+¢). Hence by (18) equation (29) for the 
asymptotic form of the secondary flow due to a source should read 


where as well as the U’s a 47 was accidentally missing from the version 
printed and a y misprinted as x. 

A still graver error is made on p. 37 in the calculation of v., the part of 
the asymptotic form of the secondary velocity field due to the trailing 
vorticity. The result given by equation (20) is doubly wrong; there should 
be a minus sign in front of this Biot-Savart integral, but more seriously the 
limiting process has been carried out incorrectly. ‘This can be seen at once 
from the fact that the derived result (22) is a velocity field which is not 
irrotational in the region upstream of the vorticity which generates it. 

The full Biot—Savart integral is 

1 1 
| dey, (82) 
and for V(1/r—r,)) in this to be replaced by V(1/r) for large r, as in (20), 
w,(r,) must tend to 0 sufficiently rapidly as r,; > «©. It is not sufficient 
that the part which is O(r7*) be removed, as was done in the paper. The 
trailing vorticity itself must be removed, since it gives a contribution which 
cannot be estimated correctly by replacing V(1/ r—r,/) by V(1/r). 

The correct answer may be obtained by subtracting a distribution of 
‘horseshoe vortices’, to remove the trailing vorticity. If v(r; v9, 2) is 
the velocity field of a horseshoe vortex with straight arms stretching from 
(~,0,0) to (0,0,0), thence to (0, yo, 29), and thence to (0, V9, 39), and unit 
circulation in the positive sense about each arm, and we put 


OX( Vo, 2 
A 20) din dey (83) 


and w, = curlvs, then w, has the same distribution of trailing vorticity as 
w,, and w,. —ws tends to zero as r > © in all directions. It is also a closed 
system of vortex lines, for the argument of equation (21) applied to w,—ws 


e. 
2S 
7 
rt 
a 
: 
r 
= 
a 
a 
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shows that [(w.—w,)d7 = 0. Hence its velocity field falls off like 7? and 
may be neglected in the present context, so that v, and v; are asymptotically 
equal. 

Now, the theory of the horseshoe vortex tells us that far from the vortex 


the curly bracket being the potential of a line of doublets of uniform stength 
(0, —3o,¥») per unit length stretching along the positive x-axis. Substitu- 
tion of (84) in (83), with use of (14), gives 


Hence the corrected form of (24) is 


AV,+V,)( v « \ AVicf y x 
(- 5.5.0) + (24) 


Of particular interest is the value of (24) on y = s = O for negative w. 
Since the curly bracket becomes y/2x? under these circumstances, we see that 
A(V,+3V, 
v,(a, 0,0) ~ 0} 
as x->— «x. The uncorrected (24) gave V,,+2V', instead of V,+3V, 
in (86). 

The above correction atfects the calculation of the asymptotic form of 
the so-called ‘tertiary flow’ (equation (25)). But this can in any case be 
criticized, as being based on a ‘ tertiary vorticity’ field derived by considering 
only the stretching and rotation of the undisturbed vorticity by the sum of 
the primary and secondary flows, and not the shearing of the secondary 
vorticity pattern by the undisturbed shear flow, which gives tertiary 
vorticity terms of the same order at large distances. Equation (26) for the 
type of equation which needs to be solved to get a better picture of the flow 
at large distances may be criticized on the same ground; the full linearized 
form of Helmholtz’s equation for the vorticity is 


(86) 


(U+Ay) = A +(Aw,,0,0), 

and the last term (representing the effect just mentioned) should not be left 
out. However, the general conclusion that there are difficulties here, 
which need to be put right by some method of the general character described, 
is right; the correct answers about the form of the disturbances at large 
distances will be presented in a forthcoming paper in this journal. 

Finally, the last formula of §3 must be recalculated in the same way as 
(85) above; it becomes 
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